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A SOLUTION OF THE GEODETIC BOUNDARY 
VALUE PROBLEM TO ORDER e 3 

by 

R. S. Mather* 

Geodynamics Branch 

SUMMARY 

A solution is obtained for the geodetic boundary value problem which de- 
fines height anomalies to ±5 cm, if the Earth were rigid. The solution takes into 
account the existence of the Earth’s topography, together with its ellipsoidal 
shape and atmosphere. 

A relation is also established between the commonly used solution of Stokes 
and a development correct to order e 3 . The data requirements call for a com- 
plete definition of gravity anomalies at the surface of the Earth and a knowledge 
of elevation characteristics at all points exterior to the geoid. hi addition, 
spherical harmonic representations must be based on geocentric rather than 
geodetic latitudes. 

No unique solution is possible in theory at the present time due to the nature 
of the Earth's atmosphere and the limited knowledge of its structure. Practical 
solutions which are only marginally in error with respect to the estimates of 

*0n leave of absence from the University of New South Wales, Sydney, Australia. 
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accuracy given above, are possible if an adequate model were adopted for the 
atmosphere. 

A quick-look analysis based on statistical considerations of the Earth’s 
gravity field, indicates that a definition which would meet the requirements 
given above for studies of sea surface topography, is afforded by a global grid 
with a 10 km spacing in non- mountainous and undisturbed regions, provided such 
information were 

(a) controlled by a global gravity standardization network of ±50 /ugal 
accuracy; and 

(b) elevations were based on a correlation of all the major continental 
datums with errors kept below ±15 cm. 

Any predictions that are necessary must be based only on the height corre- 
lation characteristics over limited distances. 
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A SOLUTION OF THE GEODETIC BOUNDARY 
VALUE PROBLEM TO ORDER e 3 

1. INTRODUCTION 
1.1 Preamble 

Until recently, it would have appeared rather inconsequential to spend time 
formulating a solution of the geodetic boundary value problem to order e 3 for a 
variety of reasons. In the first instance, others have published developments 
with this end in view (e.g., Zagrebin 1952; Molodenskii et al 1962, p. 53 et seq.; 
Bjerhammar 1962). Secondly, it seemed highly unlikely that such determinations 
could ever be put to any practical use. Further, the development of laser track- 
ing systems which promise ranges to objects in near Earth orbit with an internal 
precision of a few cm, tends to obviate any reason for carrying out the burden- 
some task implicit in the very accurate solution of the boundary value problem, 
on the basis of geodetic considerations on continents alone. The role of such 
solutions is the definition of ellipsoidal elevations, through the height anomaly, 
and hence geocentric position to a few cm. Accuracy of this type is called for 
only when studying secular variations in geodetic position which, at the time of 
writing, should be more conveniently obtained either from the ranges to satellites 
from a truly global network of tracking stations when adequate systems are op- 
erational, or from Very Long Baseline Interferometry (VLBI). 

Interest in this problem has been revived by two recent developments. Firstly, 
the definition of sea surface topography to optimum levels for oceanographic 
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analysis requires a solution of the geodetic boundary value problem to two 
orders of magnitude better than that afforded by Stokes’ integral (e.g., Heiskanen 
& Moritz 1967, p. 94) alone or one order of magnitude better than those solutions 
which took into account the effect of the topography (e.g., Molodenskii et al 1962, 
p. 118; Moritz 1966; Mather 1971b). Secondly recent developments in metrology 
promise that greater precision may well be achieved in the definition of the 
Earth’s gravity field, enabling the establishment of a global gravity standardiza- 
tion network with an absolute accuracy which is an order better than is possible 
at present. 

These proposed investigations of sea surface topography also have great 
geodetic significance in view of the commonplace departures of "Mean Sea 
Level" from an equipotential surface, as defined from the results of geodetic 
levelling, among other factors. The magnitudes of the stationary departures of 
sea surface topography, as measured at coastlines, from an equipotential sur- 
face, appear to be as large as 2 m along the north-east coastline of Australia 
(Roelse et al 1971), while discrepancies of a lesser though nevertheless signifi- 
cant magnitude, have been reported in the United States (Sturges 1972). As will 
be shown in section 4, a preliminary definition of the sea surface topography must 
precede the evaluation of geopotential differences with respect to the geoid, at 
points on the surface of the Earth, if an accurate solution of the boundary value 
problem is to be obtained, free from serious systematic error. 
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The goal of programs for the mapping of the sea surface topography from 
space seek the resolution of those characteristics with wave lengths of 200 km 
to ±10 cm (Weiffenbach 1972). This, in turn, calls for the definition of the 
equipotential surface corresponding to "Mean Sea Level" to this same order of 
accuracy over the oceans. The following development, along with all other 
means for assessing the problem, indicates that the solution of the geodetic 
boundary value problem is the most promising method available for tackling 
this problem with the accuracy quoted in the title of this paper on the basis of 
the technology available at the present time. 

While several second order solutions are available, most of these efforts 
have concentrated on amending the reference surface from a sphere to an 
ellipsoid of revolution and defining the relevant correction terms. None of the 
solutions consider the effect of the Earth’s atmosphere. Also neglected are cer- 
tain marginal conditions in the inter-relationship between the gravity anomaly 
and the disturbing potential which are of significance in defining the height 
anomaly to ±5 cm (i.e., o {e 3 h d }). The equivalent precision required in the 
definition of the gravity anomaly can be seen from equations 13 to be ±50 /j.gal. 
This figure is about four times smaller than the absolute accuracy of any of the 
stations included in the International Gravity Standardization Net 1971 (Morelli 
et al. 1971). This however does not imply that the individual values defining the 
gravity field have to be established with this precision when solving the boundary 
value problem by quadratures. This is discussed in greater detail in section 4.6. 
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A further important consideration is the preservation of geocentric char- 
acteristics of the gravitational solution. If the solution is not referred directly 
to the geocenter (Earth's center of mass), it must nevertheless be possible to 
relate the origin of the coordinate system used, to the geocenter without ambi- 
guity and to the desired accuracy. 

The development presented in the following sections, endeavors to define a 
solution of the geodetic boundary value problem with a resolution of ±5 cm in 
the height anomaly, taking into account, the effect of the atmosphere and, at the 
same time, using spherical harmonic expansions only when the function con- 
cerned satisfies Laplace's equation to the requisite precision. To emphasize 
the point, spherical harmonic functions are not used as a convenient three 
dimensional representation, but only when physically justified. Any exceptions 
to this rule are carefully qualified. In addition, the development is biased 
towards relating solutions obtained by the use of Stokes' integral alone, to that 
which is correct to o { e 3 } . This would imply that only correction terms need 
be evaluated to completely define the solution if sufficient precision were 
maintained in the calculation of Stokes' integral. These terms are formulated 
on the assumption that the solution is iterative, requiring the Stokesian term 
computation as a pre-requisite for evaluation. This procedure seems difficult 
to avoid in any solution with pretensions to accuracy, except at the expense of 
loss of definition in the context of Earth space. Section C of the Appendix shows 
the equation which needs to be solved if an iterative process is to be avoided. 
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The development also investigates techniques for optimizing the solution of 
the integral and the representations required for the global gravity field. 

1.2 Notation 

The symbols adopted have been designed to minimize confusion. To achieve 
this end, subscripts have been used to differentiate between quantities which 
have similar characteristics. Thus the symbol V is used to represent a potential 
whose magnitude is usually small. V d is the disturbing potential, while V is 
the potential of the atmosphere. Similarly the symbol h is used to represent 
ellipsoidal elevation, while h d is the height anomaly. The subscript d also 
traces a common thread, namely, quantities which are a consequence of the 
distortion of the Earth from an ellipsoidal reference model. 


1.2.1 Symbols 


A = constant associated with azimuth 

A = surface harmonic of degree n in the spherical harmonic representation 

n 

of disturbing potential 
a = equatorial radius of reference ellipsoid 
C nm = surface harmonic = P nm (sin <f> c ) [Ci nm cos mA. + C 2nm sin mXj 


(i-sin 2 ^) + ^ + o« 2 > 

' m 


+ + c„ + O 


c - 


_ cos (1/2 \p - 6) 1 

cos (1/2 \p + 6 + S) ' 


1 9 dR 

1+2 


C A = 


R 

1 + c— 


- 1 


(A-6) 
(A- 14) 
(A-27) 

(65) 
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c^f + m-3fsin 2 0 c (A . 37) 

dR = h max - h + f s in 2 <fc c + o { f 2 R} (56) 

dS = element of surface area at the physical surface of the Earth 
dS' = dS cos /3 = R 2 dcr 
dV - element of volume 
do- = element of solid angle 

E {A g} n = error of representation of gravity anomalies for a n° x n° square 
e = eccentricity of the meridian ellipse = 2f - f 2 
F (0) = f(i p) sin 0 

f = flattening of the meridian ellipse 
f ( 0 ) = Stokes* function = cosec 1/20 + 1 - 5 cos 0 - 6 sin 1/2 0 - 

3 cos 0 log [sin 1/2 0 (1 + sin 1/2 0)] (82) 

G n = n-th degree surface harmonic in the representation of A g * at the 
surface of the Earth 

g = observed gravity at the surface of the Earth 
h = ellipsoidal elevation 
h d = height anomaly 
h n = normal height 
h' = orthometric height 

K — constant for evaluation of Stokes* integral by quadratures 
= 1.58 x 10 2 cm mgal” 1 (degrees) 2 
k = gravitational constant 
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M = mass of the Earth 

M{X} = global mean value of X 

2 

a co 

m = 

7 e 

N = elevation of geoid above ellipsoid 

r = distance from geocenter to a point at the Earth's surface 
R = radius of minimum geocentered sphere which encloses the solid earth 
r = mean radius of the Earth 

m 

r = distance from the element of surface area dS to the point of computa- 
tion P at the Earth’s surface 


r = 2R sin 1/2 0 (61) 

r 0 =2R m sin 1/2 0 (A-12) 


U - spheropotential due to the reference system 

U Q = U on the surface of the reference ellipsoid, which is defined as an 
equipotential surface 

V a = potential due to the atmosphere 

V d = disturbing potential 

K = v d - V a 

W = geopotential 

W Q = potential of the geoid 

X* = geocentric rectangular Cartesian coordinate system X l X 2 X 3 

= local rectangular Cartesian coordinate system x t x 2 x 3 with x 3 axis 
along local normal, the x x x 2 plane defining the local horizon and 
completing the local Laplacian triad 


7 



11 


a = azimuth 
fS = ground slope 

/3 = term of order f in the formula 

y 0 =y e [1 + sin2 + ^2 sin4 ^ 

for normal gravity 
y = normal gravity 

Ag = gravity anomaly at the surface of the Earth 


W 0 -U 0 c)Ag 

Ag' = Ag - 2 — - dR 8 


R 


3R 


(83) 


Ag = Ag + 2-^ [f + m - 3f sin 2 0 c ] + 2 ~ \ eC 2 + dR + o{f2Ag} 


= Ag + o {1 Cf 2 Ag} 

h - h 

AR = R m [c Rp - c R ] = f (s in 2 4> c - sin 2 0 cp ) + -^ + o{f } (A- 11) 

m 

AW = difference in geopotential between the geoid and a point at the Earth's 
surface 


S = f sin 24> c cos a CT + o{f 2 } (A-10) 

t) = prime vertical component of the deflection of the vertical ( — ,c , 2 ) 

0 cot ii b - S + o{f 2 } for \p > 10° (A-26) 

2R 2 y 

m 

A = longitude, positive east 

g = meridian component of the deflection of the vertical ( = ^ ) 

£ = deflection of the vertical, positive if the vertical lies north, east of 
the outward normal 
P = density 

2R (1 + C } 

$ = _J? — [(1 + c R ) cos 8 - (1 + c Rp ) cos (0 + 8)1 - 1 (A-17) 

r 0 c 1 + C r) 
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i> c = geocentric latitude, positive north 

0 = angle between geocentric radii to the element of surface area dS 
and the point of computation P 
co = angular velocity of rotation of the Earth 


1.2.2 Conventions 

a = b + o{b 2 } = terms whose order of magnitude are equivalent to or less than 
b 2 are neglected (b < 1) 


X aVa = x iyi + x 2 y 2 

— d J y = x dy + — d 2 y + — d 3 y + 
i! 2' J .V 


i taking all possible 


values 

x i = a ij b j " x i = a ii b i + a i 2 b 2 +..... , there being as many equations 
as possible values of i 

a fi c a has the same order of magnitude as c 
a = c a is approximately equal to c 


1.2.3 Subscripts 

a = assumed values, usually either astronomical or with reference 
to a regional geodetic datum 
c = geocentric; correction to free air term 
d = disturbance value between physical and reference systems 
e = equatorial value 

g = geodetic values referred to the geocentric ellipsoid 
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m = global mean value 

p = evaluated at the point of computation P 
a - evaluated at the element of surface area dS 
X = value of "x" on the surface of the minimum geocentered sphere 
which encloses the Solid Earth. 


2. BASIC DEFINITIONS 
2.1 Gravitational Potential 

Gravitational potential is defined as the scalar W such that the acceleration 
vector g due to the gravitational field is defined by the relation 


g = - vw 

where 



( 1 ) 

( 2 ) 


the X i axis system being a geocentric Cartesian frame whose Earth space lo- 
cation is defined by the unit vectors i along the X i axes. 

Along the equipotential surface W = Constant, 



(3) 


if s is a linear displacement on the surface. If the latter is defined by the vector 
S, given by 


R = X i T, 
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equation 3 can be expressed as 



ds c)X. ds ds 


l 


(4) 


As the vector dR, as shown in figure 1, lies entirely in the equipotential 
surface, it follows that the vector VW which equals -g from equation 1, is normal 
to the equipotential surface. 

The significant conclusion is that the vector g is always normal to the re- 
lated equipotential surface. The incremental normal displacement is called an 
increment in ortho metric elevation. 


2.2 The disturbing potential 

The disturbing potential V dp at a point P in Earth space is defined by the 
relation 

V . = W - U (5) 

d P P P ' 

where W is the geopotential due to the rotating Earth and its atmosphere, and U 
is the spheropotential due to the system of reference which arises from a gravi- 
tating ellipsoid of revolution rotating with the same angular velocity as the Earth. 
This definition implies a rigid Earth and deviations from this model are dis- 
cussed in section 2.4. The subscript p refers to evaluation at the point P. 

The definition of V d at the Earth’s surface is not achieved in circumstances 
identical with those at satellite altitudes. In the latter case, V d is determined 
directly from observations. In such a case, the position of the point P at which 
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V d has been defined is therefore known. This is not so at the Earth’s surface 
where the situation is more accurately described in figure 2. The observed 
quantities specifying the Earth space location of a point P at the Earth’s surface 
are 

(a) the assumed latitude 4> a and longitude \ ; and 

(b) the difference in geopotential AW between the geoid and P as determined 
by geodetic levelling. 

and \ can be either observed astronomically or else defined with 
respect to some regional geodetic datum. The telluroid has been defined as the 
locus of points Q(0 a , k a , U Q + AW ) on the reference system where the first two 
coordinates are astronomical values, U 0 being the potential on the surface of the 
reference ellipsoid (Mather 1968, p. 518). hi the context of the geodesy of the 
70 ’s, it is more likely that <p a and k a are coordinates on the regional geodetic 
system, which differ from the equivalent values (0 , X. ) on a geocentric system 
by up to 5 arcsec. More about this is section 2.3. The elevation h of P above 
the ellipsoid is not known, but that of Q is. U p which is therefore unknown, can 
be related to the value U Q at both Q and P', situated at the intersection of the 
normal through P and the equipotential surface U = U Q , by the Taylor series 



BMJ 

Bh 1 


(6) 


where h d is the height anomaly at P (= PP’), measured along the spherop 
normal at P. The latter deviates from the ellipsoid normal by an angle whose 


12 



16 


magnitude is of order f 2 (Mather 1971b, p. 80). It follows that the effect of 
curvature of the normals introduces linear errors into elevation whose maximum 
order of magnitude is f 4 x 10 km. The linear equivalent is 10 2 mm and of no 
consequence in the present development. 

As h d is therefore normal to the spherop U = U Q , it follows from equation 1 

that 


BU 

Bh 


= - 7 


(7) 


where y is normal gravity at P'(0 g , , U 0 + AW ). It should be noted that 


7 = 7 q - (</> a " 0 B ) T e sin 2<£ c 

= 7 q + A£ y e ySj sin 2<£ c + 0(1 Mgal} 


( 8 ) 


y e and /3 1 being the relevant terms in the formula for normal gravity (e.g., 
Heiskanen & Moritz 1967, p. 78), while A^f is the correction to the meridian 
component of the deflection of the vertical due to the departure of the regional 
geodetic datum from a geocentric location (e.g., Mather 1971a, p. 63). Equation 
8 is of relevance only if a world geodetic system is not available. As the defini- 


tion of y is only required to ±50 7 gal, it would suffice if A^ were resolved to 
approximately ±1 arcsec (as a prerequisite to a complete solution) in such a 


case. 
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The second derivative is well known to be 


<5h 2 


= 0.3 mgal m -1 . 
oh 


Thus the term obtained in equation 6 when i = 2 has a magnitude of 
o{10~ 3 kgal m} and can be neglected in the present study. Equation 6 can 
therefore be expressed as 


U p = U Q + AW - yh d + o{10“ 3 kgal m} (9) 

The geopotential W p at P is unknown because the potential W Q of the geoid 
has not been established. Also 


Wp = w 0 + AW (10) 

The use of equations 5, 9 and 10 give 

V dp = W o - u 0 + ^ h d + 0 { 10 " 3 k g al m} (11) 

In summary, the height anomaly h rf is the linear displacement along the 
spherop normal, of the geop W = W p passing through P at the Earth's surface, 
from the associated spherop U = U Q which has the same difference of potential 
with respect to the reference ellipsoid U = U as W = W has in relation to the 

0 p 

geoid W = W . Thus 

o 

U q- U o= W p- W o= AW - 
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If W = U„ , W = U . This cannot however be assumed to be the case at this 

0 O p Q 

stage. 

Notes 

(i) y is the value of normal gravity at a point on the associated spherop 
U = U Q . Its relationship to the value of normal gravity at an equivalent 
point on the reference ellipsoid is defined by equation 96, the required 
precision being o{ e 3 y] . 

(ii) The term W Q - U 0 is indeterminate from gravitational considerations 
alone. It can be evaluated if a geometrical relation is established 
independently between the geoid and the ellipsoid. Its magnitude has 
been estimated at 2.7 kgal m (Mather 1971b, p. 98). Else it can be 
assumed to be zero on the basis that kM has been determined to 
o{4 x 10 12 cm 3 sec -2 }. Present day determinations (e.g., Esposito 
1972) claim an accuracy of 500 x 10 12 cm 3 sec -2 and hence fall short 
of the precision required for satisfying this condition at the time of 
writing. 

2.3 The gravity anomaly 

The gravity anomaly at the surface of the Earth Ag is defined as the differ- 
ence between observed gravity g p at P on the Earth’s surface, and situated on 
the spherop W = W p , and normal gravity y at the point P' on the associated 
spherop U = U Q , as shown in figure 2. Thus equation 8 gives 
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Ag = g p - y = g p - r Q - r e Pi s in 2< ^c 

or 

Ag = Ag a - y e sin 2 4> c ( 12 ) 

where Ag a is the gravity anomaly calculated using geodetic coordinates referred 
to the local geodetic datum. A£ is given by (Mather 1971a, p. 63) 


Acf ~ 1 ■ [A£ 0 (p 0 + h ) [cos 4> 0 cos 4> + sin sin^cos 8\] 

p + h 


+ At) 0 (v q + h Q ) sin 4 > sin 8 A - 


- AN Q [sin 0 O cos 0 - cos 4> 0 sin <f> cos S X.] 


where p , v are the radii of curvature of the reference ellipsoid in the meridian 
and prime vertical directions, while A^ 0 , At^ 0 and AN 0 are the corrections to the 
deflections of the vertical and the ellipsoid elevation at the origin of the regional 
geodetic datum, on conversion to geocentric values, the subscript 0 referring 
to values at the origin, and S \ = k - A 0 . As V d is defined as 


V d = w - u, 

differentiation along the spherop normal gives 

_ BW BU 

Bh Bh Bh ’ 
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all values referring to P at the Earth’s surface. From equations 1 and 7, 

clU 

"3h = r p’ 

while observed gravity 

aw 

g ~ “ ah' ’ 

where differentiation in this case is along the local vertical. Small changes in 
h’ are quantities which can be observed, but h’ itself, which is the orthometric 
elevation, is unknown in the absence of knowledge of the stratification of matter 
exterior to the geoid. Thus 


Bv d 

— = - g cos £ + 7 p > 


where £ is the deflection of the vertical. y p is not a known quantity while y , 
as defined in equation 7, is. y p can be related to y by a Taylor’s series 


h d 

y p = r + Ty 


Bh 1 


a 2 y/ah 2 is given by (e.g., Heiskanen & Moritz 1967, p. 79) 


— — = — = o{2 x 10” 14 cm -1 sec -2 } . 

ah 2 " a 2 


The term obtained when i = 2 will have a maximum magnitude for the largest 
possible value of h d (= o {10 2 m}), given by 
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A- h d ^-2. - o {1 /xga 1} . 

2 Bh 2 

Thus 

r— = - g 1 - — £ 2 +y+h, ^2 + o {1 ugal} 

Bh L 2 J d Bh 

= ~ A S + h d g i 2 + o {1 /xgal} (13) 

as C = 0 {3 x 10~ 4 rad} in mountainous country. 

Notes 

(i) In estimating magnitudes of quantities, Ag should be assigned 

o {10 2 mgal}. Thus e 3 Ag = o (5 x 10/j.gal}. The contribution of the 
term 1/2 g£ 2 holds the same sign at all locations with a maximum 
magnitude of o {5 x 10 /^gal} and must therefore be treated as a 
systematic effect. It will be retained in all formulae for the present. 

2.4 The Boundary Value Condition 

The formulation of the boundary value condition is freely available in the 
literature (e.g., Moritz 1965; Mather 1968). Derivations stem from Green's 
third identity (e.g., Heiskanen & Moritz 1967, p. 11). If r is the distance of the 
relevant element of surface area dS or volume dV. interior to the bounding 
surface S, from a point P on S, the scalar 0 satisfies the equation 

f i V 2 0 dVj .= - 27 T<t> p + jjJ* jl V • N<£ - <£V • N I j dS (14) 
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where N is the unit vector defining the outward normal to S, V being defined by 
equation 2. No approximations are involved in the formulation of equation 14 
apart from assumptions implicit in qualifying the existence of the relevant inte- 
grals. v • is the derivative of the scalar as evaluated along the outward 
normal and must exist exterior to and on the surface. The geopotential W is 
given by 

w = v e + V a + v r (15) 


where V is the attractive potential due to the solid Earth and oceans, hereafter 

e 

referred to as that of the solid Earth, V is the attractive potential due to the 
atmosphere, and V r is the rotational potential. As V e satisfies Laplace's equation 
at all points exterior to the physical surface of the Earth S, while V does like- 
wise at all points in V. within S. To simplify the application of equation 14 to V, 
with ^ exterior to S gives 


277 V 


ep 


-IE 


1 V • N V - V V-N — dS 

e e f 


(16) 


Similar application to V r with interior to S gives 


2 a? 


f IdV. + 277 V rp = J^* [iv-NV r -V r V-fii] dS 
V . s 


(17) 


where co is the angular velocity of rotation of the Earth, which is assumed con- 
stant, implying a rigid Earth. It also follows that 


19 



23 



Similar application of equation 14 to the spheropotential U due to the gravitating 
reference system which has the same rotation characteristics as the Earth, 
and defined by 


U = U + V 

e r 


(19) 


gives 


2w ( U ep- V rp) = -JJJ [^*N(U e +V r )-(U e + V r )V-Nl 


dS 


+ 2 or 


!>■ 


( 20 ) 


As the integrations in both equations 18 and 20 are taken over the same sur- 
face, it follows that appropriate differencing gives 


2 -(W p - U p - V a ) 


= -J[ [i 9 'N(V c -U.)-(V. dS (21) 


on using equations 15 and 19. Further combination with equation 5 gives 


27f ( V d P - V a P > = - jj [j^NV* -V'V-nI] dS (22) 
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where 

v; = v d - V, (23) 

Notes 

(i) V' is the disturbing potential due to the solid Earth and oceans. It is 
of importance as it satisfies Laplace's equation at all points exterior 
to the Earth’s surface provided that the ellipsoid lies within the former 
at all points (i.e., the ellipsoid is smaller than the geoid by at least the 
maximum negative geoid undulation). This would require in theory, an 
ellipsoid which is approximately 100 m smaller than that of best fit. 
Under these conditions 

V 2 v^ .= 0, 

and hence V' can be represented by a solution in spherical harmonics 

d 

of the form 

oo n 

df E c - < 24 > 

n“0 m = 0 

where 

C = P (s in 6) [Ci cos mk + C 2 sin m\] (25) 

(ii) The harmonic of degree one can be excluded from equation 24. This 
would imply that the reference ellipsoid were centered at the center of 
mass of the solid Earth which will not coincide with the geocenter unless 
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the mass distribution of the atmosphere has no first degree harmonic 
on a geocentric coordinate system. Further, no unique solution of the 
boundary value problem is possible unless the density distribution of 
the atmosphere were known. As this is constantly varying, the definition 
of a model for the atmosphere is called for. 

Let G e and G a be the centers of mass of the solid Earth and 
atmosphere respectively, as shown in figure 3. If the X. axis system 
is centered on the geocenter G, it is possible to define the coordinates 
X. of G e in terms of those (X ai ) of G a , on the assumption that the 
model for the atmosphere is capable of formulation from direct 
measurement more readily than that of the solid Earth. 

If dV is an element of the volume V e exterior to the solid Earth 
which contains the atmosphere with density p , then 

M . X., = JjJ p. x i dv , < 26 > 

where M a is the total mass of the atmosphere. A similar consideration 
of the solid Earth, of total mass M e , contained within the volume V. , 
gives 

M , X „ = JJJ P. x i ^ (27) 

P e being the density of matter contained in the element of volume dV A . 
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As the X i axis system has its origin at the geocenter, 


Mx 


0 = J J J p. x i d V, + J JJ P' x, d V, = M, X„ 


+ «e 


where M is the mass of the Earth. Thus 




(28) 


In summary, the boundary condition set out in equation 22, is built 
around the disturbing potential for the solid Earth (V') which has the 
advantage of satisfying Laplace’s equation and hence being expressible 
in spherical harmonics. This representation would not have any terms 
of degree 1 if the atmospheric potential, referred to a geocentric co- 
ordinate system, also had no first degree terms. If this is not the case, 
as seems likely, the spherical harmonic representation is referred to 
a coordinate system based on the center of mass of the solid Earth, the 
relationship to the geocenter being given by equations 26 and 28. This 
problem will not be considered further. 

(iii) The inclusion of equation 17 in equations 18 and 20 implies that eo is a 
constant independent of position within the surface of the Earth. This 
is not so in practice due to departures of the Earth from a rigid body, 
variations in the rate of rotation and polar motion. The first effect is 
allowed for as the well known correction for Earth tides to observed 
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gravity, with a magnitude o {lOV gal}. The gravitational effect of the 
rotation g r is given by 


g r = POJ 2 


where . 


P = 


(X? + X*) 


1/2 


the Xj being defined as in section 1.2.1. The change Sg r in g r due to 
polar motion can be interpreted as a consequence of changes dXj in 
X i and d co in co, the relevant relation being 


X. d X, + X 2 d X. 

S g r = co 2 + 2 p co d co 


X a dX a dw 

+ 2 

L p 2 " 




As d co/co = o (3 x 10 7 } and dX/X = o (10 6 } at mid latitudes, the ef- 
fect on g = o { lp. gal } . The limited magnitude makes it possible for 
this effect to be neglected for the present development, even though it 
is do minat ed by a set of even zonal harmonics. 

(iv) The validity for adopting spherical harmonic representations for func 
tions on S is discussed in section B3 of the Appendix. 
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3. SOLUTION OF THE BOUNDARY EQUATION 
Equation 22 can be written as 



the surface integral being taken over the physical surface of the Earth. The 
latter can be represented by the telluroid without introducing errors in excess 
of o{ f 2 } and hence smaller than the accuracy sought in the present study. On 
adopting a local x . axis system at the element dS, with the x 3 axis oriented along 
the local spherop normal and the x x and x 2 axes oriented north and east respec- 
tively, it can be shown (Mather 1971b, p. 80) that 
/ 

N = cos /3 [- tan f3 a a + 3] (30) 

where /3 is the slope of the topography at dS while / 3 1 and fi 2 are the components 
of the ground slope in the north and east directions respectively. Thus 

V-N 1 = cos £ [x a tan /3 a - x 3 ] (31) 

r r 3 

and 

V-NVj = cos fi 

On considering equations 11 and 23, 


d x 


tan 


t + 3x. 


(32) 


V d = V d - V a = (W 0 - U 0 ) - V a + 7 h d + o (5 x IQ' 2 kgal m} (33) 
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The use of equations 13 and 23 gives 

3 Vj 3 Vj 3 Vj _ 3 v - 
3 h 3h 

= -Ag + h d ^Z+Ig^ 2 - _1 + o {1 fj. gal} 


(34) 


The last term in equation 34 is the attraction of the atmosphere. It follows 
that observed gravity must be numerically increased for the attraction of the 
atmosphere before computing the gravity anomaly. Assuming a density of 10" 3 
gm cm -3 for a 20 km thick layer, treated as a Bouguer plate, the magnitude of 
this term is 


3 v « 

— — = o (10 3 gal} . 
dh 

If standard concepts of the nature of the lower atmosphere are accepted, (e.g., 
Smithsonian Meteorological Tables 1958, p. 267), this correction will be corre- 
lated with elevation and is more than likely to have a first degree harmonic as 
discussed earlier in note (ii) to section 2.4. For a treatment of the atmosphere 
consisting of a series of nearly ellipsoidal shells, see IAG 1970, p. 62 et seq. 
The use of equations 31, 32 and 34 in equation 29 gives 


v dp = V. 


ap 




cos/3 dS 


(35) 
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Cos /3 dS is the projection of the element of surface area dS from the 
telluroid, onto the associated spherop (Mather 1971b, p. 81). The projected 
surface area is related to the element of solid angle do- by the relations 

d S cos /3 - R 2 d a = R 2 cos <£ c d 4> c d k ( 36 ) 

where (R, 4 > c , are coordinates on a geocentric spherical system. 

The basic equation at 35 can be written without approximation as 


where 


and 


= v„„ + I + L 


dp 


ap 


(37) 



The integral I A contains the standard Stokesian term, but masked by 

(a) the ellipticity of the meridians; 

(b) the undulations of the topography; and 

(c) the gravitational effect of the Earth's atmosphere. 

1/r can be expressed as the standard zonal harmonic series (e.g., Jeffreys 
and Jeffreys 1962, p. 634) 


1 

r 


n — fl 


( 40 ) 
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As Vj satisfies Laplace's equation at all points exterior to and, in the limit, on 
the Earth's surface, as discussed in section B3 of the Appendix, can be ex- 
pressed by the series 


= n 

L—j R" + 1 


* 1 


(41) 




The exclusion of the first degree harmonic places the reference ellipsoid at the 
center of mass of the solid Earth, which does not coincide with the geocenter. 
The resulting consequences have been explained in section 2.4. It follows that 


3 V d 3 V ' 

— = — — sec (S 0 - 8 0 ) 
Bh BR v c ' 


uu 

-L 


(n + 1) 




+ o 


n = 0 


>n + 2 



(42) 


where 8 0 and are defined in figure Al. The expression for the gravity 
anomaly A g in terms of spherical harmonics is obtained from equation 34, on 
defining the vertical gradient of normal gravity at the surface of the Earth. This 
can be related to the equivalent value ((3y/Bh) 0 ) at the ellipsoid by the Taylor 
series 


B y _ 
Bh “ 



The first term on the right is well known to be (e.g., Heiskanen and Moritz 1967, 
p. 293) 
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(lZ^ = .111 (1 + f + m - 2 f sin 2 <P c + o {f 2 }), 

\3 h / a 

0 

where all quantities are as defined in section 1.2.1. As y Q on the reference 
ellipsoid is related to 7 at the surface of the Earth by the relation 

7 0 = 7 ^l + 2-|j-+o {f 2 } j , 

and a is related to R through equations A3 and A4 as 

a = Rq (1 + f sin 2 0 c + 0 {f 2 }) = R ^1 - — + ^ s in 2 ^ + ° - 

it follows that 


= (l + f +m-3— - 2 f sin 2 <£ ^ = - — (1 + f + m - 3 f s in 2 $ + o {f 2 } ) ( 43 ) 

3 h a \ a / R 

as B 2 7 / 3 h 2 is given in section 2.3 as 67 /a 2 . 

The combination of equation 43 with equations 33 and 34 gives 


a v„ 


d V o 1 u v 

' B = *3TH (V ^ + V *- (VU » ))(1 + f + m ' 3 1 sin2 ^ ) + 5 gJ2 ‘Tif 

( 3 V # 2 V # \ 2 V' o 

- (f +m-3f sin 2 0 c ) + |(W o -U o )(l + f + m-3f sin 2 ^) 


[2 V o V \ , 

-hr + nrH BC2 + °< f!AB} 


(44) 
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The contribution of the second set of terms is o {100/xgal} while that of the 
third is o {1 mgal} if (W 0 - U 0 ) is o {10 kgal m}. The terms of o{f } in this set 
can therefore be disregarded without introducing errors in excess of o(10Vgal} . 
The use of equations 41 and 42 gives 


A g =£<"-D A " ^ 


n = 0 


Rn + 2 R 


c* + 2 


W o‘ U o 2V 3 V 


R 


R 3 h 


+ g £ 2 + o (f 2 A g} , n / 1 


(45) 


where c^ is given by equation A3 7. 

On using equations A15, 41 and 42, the integral I A at 38 can be written as 


*A " I A1 + J A2 


(46) 


where 


and 


\n = 0 K n =»0 K 

1 f fR 2 2 n + 1 A n 

n = 0 K 


dcr, n j. 1 


A2 


-4J'r 


d cr 


(47) 


(48) 


The second integral deals with quantities which are due either the topo- 
graphy or the ellipticity of the meridian. The contribution from this integral is 
f times smaller than that from I A1 for the same region, except when is small 
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and in those instances when <hp - h) is large in comparison to r. Even in regions 
of rugged topography situated very close to the point of computation P, the con- 
tribution of I A2 is at least an order smaller than that of I A1 as can be seen from 
the discussion linking equations A17, A18 and A19. Equation 47 contains the 
major contribution to I, including that due to the well known integral of Stokes. 
The spherical harmonic model adopted for in equation 41 is a necessary 
intermediate in the combination of the effects of the gravity anomaly and the 
disturbing potential. This provides an effective technique for obtaining a first 
approximation for V ' through Stokes’ integral with adequate accuracy, hence 
permitting the use of an iterative method of solution of the geodetic boundary 
value problem. For further discussion on the possibility of using non-iterative 
procedures, see section C of the Appendix. 

The conventional procedure due to Stokes cannot be followed when solving 
equation 47 without introducing approximations due to the following reasons. 

1. R varies with da. 

2. The spherical harmonic expansion 


CD 

L 

n = 0 


2 n + 1 


R n+ 2 


n * 1 


only holds at and exterior to the surface of the Earth and is defined at 
all points on the latter. In this case, the values at the Earth's surface 
are defined for limited ranges of R given by R = R ^ + o {f }such that points 
on it are completely defined by the set {4> c , k}. 
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The standard Stokesian practice calls for the replacement of the spherical 
harmonic term by a set of surface harmonics G n , implying the relation 


w CD 


(n - !) A n 


n= 0 


n ■ 0 


R 


n + 2 


n ^ 1 


(49) 


This technique which is valid on the surface of a sphere, next equates individual 
harmonics of degree n by the relation 


A 

n 


R n * 2 
n - 1 




(50) 


The method breaks down at the physical surface of the Earth as the variations 
of R from R m , though of o {f }, are nevertheless functions of the set {<p c ,k} . In 
this case, A g’, defined at equation 49, is related to the gravity anomaly through 
equation 45 as 




(n- 1) 


R 


n + 2 


2 V' 

n ^ 1 = A g + — - c Aj 


(51) 


being given by equation A37 and ca b by equation A44. Note that Ag’ is not 
defined for all {4> c , X.}, given R, unless R > R, where R is the radius of the mini- 
mum sphere which is exterior to the solid Earth, with its center at the center of 
mass of the latter. 

It can be concluded that the replacement defined in equation 50 is valid if 
the surface harmonic expansion of Ag T (a) refers to the sphere of radius R, and 
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not the physical surface of the Earth; and (b) does not exist in those regions 
where R is less than the geocentric radius to the local topography. 

Any other interpretation would result in a loss of definition, causing approxi 
mation errors of order fh d , which' is unacceptable in terms of the accuracy esti- 
mates specified for this paper. As can be expressed by equation 41 at all 
points exterior to the physical surface of the Earth, it follows that Ag* also 
exists everywhere in this same region, being defined by the first equality at 51. 
It also has the characteristic of taking values at the physical surface of the 
Earth defined by the second equality at 51 . These values can be deduced from 
observations at this surface, where A g’ differs from the gravity anomaly A g 
by magnitudes which are of order 1 mgal. 

Equation 47 can therefore be written as 


T=_l f T*I T 2 n_jl G dcr 
A1 2 tt J J r L-i 2 (n - 1) n 

n = 0 


(52) 


without approximation. It should also be noted that 



(53) 


again without approximation. The Stokesian procedure calls for the expansion 
of l/r in a set of zonal harmonics, using equation 40 and the use of the ortho- 
gonal properties of surface harmonics when I can be transformed to 
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l M = 


1 ff y 2n+l K\" 
4 77 J J n - 1 \ R / 


(cos i/i)ig' do- 


This manipulation is not possible in the case of equation 52, without introducing 
errors of order fl A i for the following reasons. 

(1) R does not remain constant as {<£ c ,A} varies over the Earth’s surface. 

(2) While surface harmonics retain orthogonal properties on integration 
over any closed surface as illustrated in section B2 of the Appendix, 
this does not apply to G n , which may not exist for certain (<£ c , iV} as 
explained in the discussion following equation 51. 

The conditions for the recovery of Stokes' integral I g , correct to o{f 2 I}, 
can be obtained as follows. As Ag’ and V ^ exist in the domain exterior to the 
physical surface of the Earth, satisfying equations 51 and 41 respectively, let 
Ag» and V ^ take values Ag 7 and V/ on R, where 


R = a + h max (55) 

h ma* bein S any number marginally larger than the maximum ellipsoidal eleva- 
tion possible. The displacement dR along the geocentric radius between a point 
at the Earth's surface and the sphere of radius R is given by 

d R = R - R = h max - h + a f s in 2 0 C + o {f 2 R} = o {f R} (56) 

V d 311(1 A g are related to VJ and A g' at the surface of the Earth by the Taylor 
series 
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Ag^Ag' + (l^>i = A g' + c 


1 ! 


(57) 


d R 1 


and 


vT = v- + iiSl ll^ = v: + c 

a d j | _. ¥ d T v 


1 ! 


(58) 


d Ri 


The use of equations 57 and 58 in 53 gives 


A1 

where 




(i + c A) |s? + !^-(c i + i^ + !^ 

2R \ 2 R / 2 R \ R 


— + o{f 2 >l 


d a 


(59) 


and 



(60) 


r = 2 R sin - di 
2 


(61) 


Equation 59 can be further partitioned according to the relations 


*Al “ I S + *c 


where 


and 




d a 


(62) 


(63) 




3 V' 


c A ,A g . + _^).( Cs + iJ.i| «♦.(*>! 


dcr+ O {f I c ) 

(64) 


On appreciating the analogy that exists between equations A12 and 61, c A can 
be shown to be 
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where 


1 + 2-^5-+ o {f 2 } 


c A=- 


R 

(1 + c T ) 


- 1 


1/2 



d R + d Rp 
R 


+ o {f 2 } 


(65) 


(66) 


on lines similar to those used in the derivations of equation A14. 

The first set of terms on the right of equation 64 contribute o {fl g } except 
when ^ is small and A R is do minat ed by the magnitude of h p - h. In such a 
case, considerations similar to those expressed in formulating equations A17 
to A19 apply. 

The second set of terms is also of order { f A g }. While this is obvious in 
the case of the last of the terms in this set, it is also apparent in the case of the 
other two. The largest contribution to c is due to the term obtained in equation 
57 when i = 1. 


c 

g 


= dR lr *°{ <dE>2 


B 2 (Ag) l 
3R 2 


(67) 


where 3Ag'/3R is well known to be given by (e.g., Heiskanen and Moritz 1967, 
p. 117) 


3Ag' 

3 R 


= - 7 



3£ a tan <£ c 2 h d f 9Ag] 

+ o -s f r 

3x a R r 2 l 9RJ 


( 68 ) 
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The magnitude of this term is o{10~9 gal cm -1 } , giving c g = o {1 mgal} if 
^ £ a /^ x a = o (l arcsee (10 2 km) -1 }. In disturbed areas, this could be one order 
of magnitude larger. In all circumstances, terms involving c g need only be 
evaluated to the order of the flattening, to meet the accuracy requirements called 
for in the present development. 

Similarly, c v is given by 


3 v d f V; 

c v = d R + ° i Cd R) 2 - 

3h 1 3h* 


(69) 


(3 V,' /dh) is given by equation 34 and the effect of the term containing c^ is f 
times smaller than the contribution of (3 V,' /3h) in the Stokesian term. Equation 
64 can therefore be written as 


Ic ^Ht ( A8> ( CA+ f tt) + ir( Cfl+ 3 ! r)- dE ^ +0<f!Ae ' } ) d 

(70) 

The solution of equation 63 is well known but will be traced out here for 
completeness. From equations 41 and 51, 1 can be written as 


O' 


_ 1 f f R 2 V 1 2 n + 1 A n , 

jjtz;— s — 


(71) 


As 


Ag ' =(A6 'VH=£](n-l)-l-=^G„, n/1 

n = A -K 

the replacement 


(72) 


n = 0 


n 25 0 


37 


41 


K = 


R n + 2 
n - 1 


(73) 


is valid and I s can be written as 




JLg 

4 77 


J r 2n pr, 

J 

0 *0 


2 sin 1 i// cos ii/i d ^ d a 

2 s in i- 1 1) 

2 


R 2 f f 2 n + 1 - 

47 J JZ_ 177ri- G “ 

n= 2 


X 


1 

R 



P. o (cos yp) da 


(74) 


on using equation 40, R p being equal to R from the definition of r at equation 61. 
As R, R p and G n exist over the range of integration at all points on the surface 
of the sphere of radius R, the use of the orthogonal properties of surface har- 
monics gives 


I = - R (M {A g'}) 


_ ,5_ ff 

i?=R + 4 7 T J J 


f (i/») A g' d a 


(75) 


where 


f 


n= 2 


2 n + 1 
n - 1 



P n0 ( COS ’/O 


(76) 


Ag» is related to the gravity anomaly Ag through equations 51, 57 and 67 as 


A g' = A g + 


2 V' 


R c 4> ~ c Ag + + ° {f2 


(77) 


= A g + o (1 mgal} in undisturbed regions 
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The Stokesian manipulation is well known to be obtained as follows 


f OA) = tj + t 


where 


(78) 


and 


'R 


t i - 2 7 (^j P n0 (cos 0) = 2 R( i-I COS 0 

“ X 7 \ R R 2 


°° i /R \ n 

t2 = 3 Z]~(lf) p «° (cos ^ ) = ^ 

n= 2 R 



'R n = 2 


R 

~j P n0 ( COS d R 

R . 


(79) 


2 r 5 fi.i.!k 

R 4 \ r R R2 


cosi/f d R 


(80) 


on using equation 40. As 


r Rdi = r r -^ 

Jr 7 Jp r 


COS 0 __ 


d R + R cos 0 

p r 


f r + R - ] 

Jr 7 (f+ R - 


- ^ cos 0 


Rp cos 0) 


d R 


-] 00 

- L r J + R cos 0 

R p ” 


J r 03 _ 

| d R + d r 

r r + R - R p cos 0 


- [r + R p cos 0 log (r + R - R cos 0)]Z, 


<=4 

R 


_ — _ r+R-R cos 0" 

C 7 — R) + R p cos 0 log _f! 

R 
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Further, 




i- 

/ 

/ \ 

2 „ \l/2 

Lim 

( r - R) = Lim 

R 

*♦( 


1 - 2 — - cos 0 ) - R 

R - 00 

R-oo 




R / 


- R cos 0, 

p 


and 


Lim 

R - 00 


log 


r + R - R^ cos i p 


R 


= log 2 . 


Thus 


R 


R - R cos 0 - r - R cos 0 log 

p r p 


7 + R - R p COS 0 
2 R 


(81) 


and f(0) is obtained from equations 78, 79 and 81 as 


1 1 


R 


f (0) = 2R [ — cos0 + 5 , 1 R-R cos 0 - r -R p cos01og 

J R R 2 / R V 

= cosecl-0 + 1 - 5 cos 0 - 6s in 1.0 - 3 cos 0 log jsin ^0 
2 2 \ 2 


r+R-R cost// 

p ' 

2R 


1 + s in -1 0 


(82) 


as R p = R and 7 is given by equation 61. As W Q - U 0 is unknown, it is preferable 
to separate it from A g ’ before computations. This is convenient as Stokes’ 
function is insensitive to zero degree effects. On defining Ag c as 

W - U 

Ag c =A7-2 -2_i (83) 

R 

where Ag* is given by equation 77, I g can be written from equation 75 as 
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I s = 2 (\ -U 0 )-RM{Ag c } + -p- f (VO Ag da 




(84) 


without introducing any approximations in the surface integral term, where 


A g = A g + 


2 V' 


~ c A g + dR + o {f 2 A g} 


(85) 


from equation 77, c^ being given by equation A3 7 while c Ag is obtained from 
equation A44 as 


1 V 3 V 

c A g = j g ^ 2 " 2 J ~ ^17 + ° < f2 A S> (86) 

Notes : 

(!) A g c is approximately equal to the gravity anomaly Ag. The second and 
third terms in equation 85 are both of order f Ag and hence do not have 
magnitudes in excess of 300 m gal. The magnitude of the final term de- 
pends on the variability of the Earth's gravitational field in the locality. 
On the average, its magnitude is of order 1 mgal, though it could be one 
order larger in regions of rapid change. As this term cannot be evalu- 
ated until the terms ( 3 «f a /Bx a ) are known, the contribution of the 
Stokesian term is therefore determined by iteration. The solution need 
be iterated only once and a convenient set of formulae for this purpose 
is given in section 4.1. 
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(2) No difficulty should be encountered in computing the atmospheric cor- 
rections whose magnitudes are controlled by the model adopted for the 
Earth's atmosphere. This correction, which approaches 1 mgal (IAG 
1970, p. 62) should be applied as part of the routine when computing the 
gravity anomaly. 

(3) It would not be adequate to use the free air reduction (e.g., 0.3086 mgal 
m -1 ) in computing the gravity anomaly. Instead, the relation 


y - y - — - 1 + f+ m- — — h - 2 f s i n 2 c£ + o { f 2 }| 

0 a 2 R n 


(87) 


should be utilized when computing normal gravity at the associated 
spherop. The quantity h n , called the normal height, is obtained from 
the observed differences in geopotential AW using the equation (e.g., 
Heiskanen and Moritz 1967, p. 171) 


h 


n 


A W 

■^T 


l + (l + f+ m- 2f sin 2 cp c ) 


AW 

~o 



+ o {f 3 } 


( 88 ) 


AW having the same significance as in section 2.4. 

(4) The first iteration for h d will be obtained from Stokes' integral as be- 
fore. This contribution, equivalent in significance to that provided by the 
free air geoid in present day solutions to the order of the flattening, 
need be calculated only once if (a) the gravity anomaly A g, computed 
to o(e 3 Ag} using equations 87 and 88, were corrected for the atmospheric 


42 



46 


effects prior to evaluation; and (b) the radius of the sphere were taken 
as R and not R . 

m 

(5) The initial iteration should also include the evaluation of the components 
of the deflection of the vertical, using the Vening Meinesz integrals 
(e.g., ibid, p. 111). The second iteration need only be the correction 
terms which are more conveniently included elsewhere, as shown in 
section 4.1, as the order of magnitude does not exceed 0.3 kgal m. 

(6) It is tempting to introduce a function of the type 


v ;o-s) = 23 

n«0 



n/1 


in an attempt to combine the effects of I A1 and I A2 to give an integral 
of the form 


j a = 4^ jpWO A e a d °-> 

where 


Ag a = Ag (1 - S'). 


This technique is not unfamiliar in the literature but is not used in the 
present development for two reasons. 

(1) There is no physical validity for defining V' (1 - S) as a spherical 
harmonic series unless (1-5) satisfies Laplace's equation. 

This cannot however be proved to be the case at the physical surface 
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of the Earth where the height anomalies are to be defined. Further 
S is a function of both the topography and the Earth's ellipticity. 

It has no definition except at the surface of the Earth. Thus while 
S can be completely defined by a set of surface harmonics, it is 
invalid to equate it to a set of spherical harmonics. The differenti- 
ation between the two cases is important as the definition of the 
gravity anomaly from the disturbing potential is implicitly based 
on the existence of radial derivatives of the latter. This follows 
only in circumstances where the spherical harmonic representation 
has a physical basis. 

(2) Difficulties are posed in interpreting the location of the center of 
the reference ellipsoid as the harmonic n = 1 is inadmissible in 
the solution. Also see section 4.2. 

4. SOLUTION SUMMARY AND DISCUSSION 
4.1 Equations for a Solution to Order e 3 

From the development in section 3, the height anomaly h d is obtained from 
equations 11, 37, 38, 48, 70 and 84 as 


h . = N f + N 

dp rp cp 


(89) 


where 


N f = 

fp 


»0- U 0 R 


- R 


y 

' n 


y 4t7 y 

' p / p 


!■ 


f 0/0 Ag r do- 


(90) 
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and 




tan ^ + 




^x tan B 

a “ a 



3 




- dR 1ST + Ag ' ( CA + 1 t) + ofe3 Ag) 


da 


(91) 


Certain terms have been adjusted in equation 91 on the understanding that 
N cp = o (l(f 1 N fp } , thereby ignoring effects of order fR in its formulation. 

The constituent terms in equations 89 to 91 are defined by the following relations, 
the equation numbers referring to their identification in the text. 


while 


v; 


= V, - V. 


Ag c = A gl + Ag 2 


where the use of equations 85 and 86 gives 


3V V 

Ag, = Ag + + 2 — 

* 3h R 

m 

and 

Ag2 = ~R^ C * + dR lif + ° {e3 Ag} 

rn 

Also 


BAg _ 


y < . ^ tan *c , N f , „ l 3a 8 \ 

L-j 9x a R R 2 l J 

a- 1 m 


(33) 


(92) 


(93) 


(94) 


(67) 
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c, = f + m- 3f s in 2 0 c + o{f 2 } 


(A-37) 


m being defined in section 1.2.1. The angle i p is obtained from 


0 = cos [sin 4> c sin 0 cp + cos 0 c cos 0 cp cos dX] 


(A- 8) 


where 


dX = X - X 


(A-7) 


The geocentric distance R is related to the mean radius of the Earth R m and 
that of the minimum sphere enclosing all topography (radius R) by the equations 


* = R m C 1 + c r> 


(A-5) 


where 


and 



(A-6) 


R = R + dR 


where 


dR = h - h + f sin 2 0 + o{f 2 R} 

max r c “ 


(56) 


Also 


r = r 0 (i + c r ) 1/2 (A-13) 

where 

r o = 2R ra sin ^ (A-12) 
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+ C Rp + C R + ° {f2} 


r = 2R sin — 0 

2 


1+2 — 

R , 

c A = 1 

(1 + c r )”* 


( ar \ 2 dR + dR p 


o{f 2 } 


AR = C Rp * C R 


The other quantities requiring definition axe 


$ = _ JR - R p cos (0 + 8)] - 1 
r 2 


dV' -p, 

— tan /3 a = - y£ a tan + N f — - tan /3 1 + o{e 2 Ag} 


a . Q R ,< N . , dh 

— tan P a = — C 1 + c x ) sm 0 — 

r 2 r 2 dr 


cos I — 0 - 9 


cos | — \p + 0 + 8 
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8 - -^5. cot-Li//-S + o{f 2 } if <// > 5° 
2R 2 


S = f s in 2cb cos a 

r c cr 


-p- = cos a^. tan /3j + sin a ff tan /? 2 


The quantity A g in equation 93 is given by the relations 


Ag = g - 7. 


where g is the value of gravity observed at the surface of the Earth, while y is 
defined by 


y = y 0 + Ay, 


where y 0 is normal gravity on the reference ellipsoid and 


a _ 2y ° hn r 3h 

y 1 + f + m - 2f sin 2 0 +o{f 2 } 

a 2a 


the normal height h n being obtained from the difference AW in geopotential be- 
tween the geoid and the point at the Earth's surface by the relation 


, AW , AW L , AW . 2 ,\ fc3\ 

h = 1 + (1 + f + m + 2f s in^ c£ ] +o{f J } 

7 0 a T 0 \ a ?o 7 
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Alternately, 


7 = 7 0 


2 AW 
a 


1 + f + m - — — 2 f s in 2 0 +o{f 2 }. 

2 ay 0 ■ 


(96) 


Notes : 

(i) A rigorous solution is obtained only if the reference ellipsoid always 
lies within the physical surface of the Earth. Such a figure is smaller 
than the figure of best fit by approximately 100 m. If the values of 
normal gravity were then based on this figure plus an independently 
determined value of kM, all gravity anomalies will be too small by 
o {2 x 10 1 mgal}. The linear effect in N fp is contained entirely in the 
first two terms of equation 90 as Stokes' integral is insensitive to effects 
of zero degree. 

(ii) Ag' is defined by equation 51. In the context of the note to equation 91, 
it would suffice if Ag' were taken to be equal to the gravity anomaly A g 
for purposes of evaluation to order e 3 h d . 

4.2 Procedure for obtaining a numerical solution 

The equations summarized in section 4.1 completely define the solution of 
the geodetic boundary value problem to the order of the cube of the eccentricity. 

The form of these equations and the discussion in section C of the Appendix 
indicate that a non-iterative approach to the solution is not possible as the 
evaluation of N c at equation 91 requires a knowledge of VJ which is obtained 
from h d using equation 33, and the components of the deflection of the vertical £ . 
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It is well known that N f contributes over 90% of the magnitude of h d (e.g. 
Mather 1971b, p. 89). This is equivalent to the free air geoid in solutions to the 
order of the flattening. In determinations to order e 3 h d , the same contribution 
is obtained by the use of Ag , defined by equation 93, in Stokes' integral, as ex- 
pressed at 90. Let the numerical value so obtained be N fl while the value de- 
duced from equation 33 for V d be V dl . 

The only other contributions with magnitude greater than fh d arise from 
the terms at A17 and A32, the former being of significance only when large 
topographical undulations occur near the point of computation (ibid, p. 86). As 

K = v d. v'J, 

the use of V dl in lieu of V' when computing these topographical corrections will 
result inavalueN cl for N c which is correct to o {10 _1 N c > (i.e., to ±lm). Let 

V d2= V <.I = v i +°U0 ' ! v;> 

The computation of Ag 2 , defined by equation 94, using V d2 and the equivalent 
values of (ibid, p. 88) and its use in equation 90 will give the balance contri- 
bution to h d from the expression for N f , the magnitude being estimated at fh d , 
though it could be as large as 10“ 2 h d in mountainous regions. If this magnitude 
is N , define 

f 2 ’ 

V d3=V d2+ lN I2 =v; + o{10- 2 v;}. 


50 



54 


The use of either V d 2 or V d3 or in lieu of V d ' when evaluating equation 91 
will result in value N e2 which is correct to o{10 -2 N c ) (i.e., to ±10 cm). Defining 


d4 


d 3 


+ — N , .= V' 

y c2 d 


o {e 


V'} 

v d / ’ 


equation 91 should be iterated a third time using VJ 4 for V d to give the final 
value of N , for N . Hence 

c 3 c 

h d = N fl + N f2 + N c 3 + ° {e3 h d } ( 97 ) 

* 

These evaluations must be completed on a global basis. No solutions of the 
geodetic boundary value problem to order e 3 h d can therefore be obtained from 
data restricted to a local region. 

Notes : 

(1) A complete solution requires the evaluation of N c to be iterated three 
times. As pointed out in section C of the Appendix, it is not possible 
to avoid the iterative procedure. Considerable economy could be effected 
if the number of iterations could be reduced by obtaining a more accurate 
estimate of V' after the first iteration. Unfortunately this cannot be 
achieved by the analysis of the orbital perturbations of near Earth 
satellites as results obtained to date indicate a lack of sensitivity to the 
topographical effects. 
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(2) As pointed out in the introduction, the principal need for an accurate 
geoid solution is in the study of sea surface topography. A resolution 
to ±10 cm can be obtained by just 2 iterations of equation 91. 

(3) It has been assumed that B 2 Ag /3h 2 has a negligible magnitude. This 
would be a reasonable assumption over oceanic areas, but may be a 
limitation in mountainous and gravitationally disturbed regions. Such 
an effect is of consequence only if it holds the same sign over consider- 
able extents as discussed in section 4.3. It would not be unreasonable 
to assume that the net effect is negligible for studies of the sea surface 
topography. 

(4) The magnitude of /3x a has been assumed to be of order ±1 arcsec 
(10 2 km) -1 , when the contribution to Ag 2 is of order 5 x 10" 1 mgal. 
This magnitude can be considered to be an average value (e.g., Mather, 
Barlow & Fryer 1971, figs. 3.2 and 3.3) though it could be a factor of 10 
greater. In such disturbed regions, which are characterized by short 
wavelengths in , both positive and negative values are equally likely. 
The overall effect is therefore small unless the disturbed regions lie 
close to the point of computation. It should also be noted that such re- 
gions invariably occur in areas of rugged topography. On the other 
hand, the Australian data referred to above indicates a significant 
number of these disturbed regions are not correlated with any topo- 
graphical feature. 
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(5) The evaluations of N fl and N f2 should be based on algorithms seeking 
a precision of 5 parts in 10 4 in the final values. 

4.3 The Representation of the Gravity Field 

It must be established that the global gravity field is capable of definition 
with adequate precision to afford the determination of h d to o{e 3 h d ). There are, 
in general, two techniques available for this purpose. The first is the determina- 
tion of gravity anomalies at the surface of the Earth by direct determinations 
of g. The second is the determination of the disturbing potential V d from the 
analysis of the orbital perturbations of near Earth satellites. Accuracies at- 
tained at the present time in the determination of g are controlled by the global 
gravity standardization network. It is expected that all gravity holdings will be 
converted in the near future to values referred to the International Gravity 
Standardization Network (IGSN 1971) whose absolute accuracy is estimated at 
±0.2 mgal (Morelli et al 1971, p. 5). This figure is a factor of 4 inferior to the 
±50 fi gal figure implicit in the formulae listed in section 4.1. Individual gravity 
ties to stations in IGSN 1971 can be made to ±0.1 mgal. This figure will be 
shown to be acceptable if the density of stations in the gravity standardization net 
is sufficiently high. 

To investigate this further, it is necessary to analyze the computational 
procedures adopted in evaluating the major contribution called N fl in section 4.2. 
For simplicity, this will be called the Stokesian contribution even though this is 
not strictly so in the case of a second order solution. On excluding the terms of 
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zero degree and adopting the system of quadratures for the evaluation of the 
Stokesian term, equation 90 can be written as 


N< f cn O 


R< cm ) X 7 T 2 

4777 X 3.24 X 10 4 


i J 


= K 


i j 


(98) 


where Ag u is the value of the gravity anomaly representing an n° x n" square, 


K = 1.58 x 10' 2 


(99) 


and fi = cos 0 c or sin 0, depending on the system of coordinates adopted. 

It is required that the errors e N in N f due to the adoption of the quadratures 
technique be kept to within the ±5 cm limit. The errors in each of the individual 
products being summed, could be of two types. The first is of an accidental 
nature, characterized by the subscript a and the second is systematic in 
character, denoted by the subscript s . It is well know that the magnitude of the 
latter per individual term in the summation, should be significantly smaller than 
the former as it holds the same sign over a considerable number of terms. 

In the case of the total accidental error e in N obtained from N sum- 

Na f t 

mations, the contribution e Nta from a product of the form 

t = k (j. f ( 0 ) Ag n 2 (100) 
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should not exceed e Na //N^. The total systematic error e Ns in N f due to con- 
tributions e„, from each of the terms at 100 bears the relation 

Nt S 

e Nts ^° {e Ns/ N t } 

if the systematic error persists with the same magnitude and sign for all N t 
terms. If N t = o {10 6 }, then e Nts = o {10" 3 e } . In practice, it is more likely 
that e Nts exhibits systematic error characteristics over some subset N' t of 
N t , behaving as an accidental error over the N£ larger subsets, where 



The evaluation of a surface integral by quadratures calls for the sub- 
division of the surface into infinitesimally small elements, the evaluation of the 
kernel of the integral at each of these elements, and the summation of each of 
these magnitudes. In evaluating equation 100, it is necessary to adopt values 
for (x , f(0) and Ag to represent the m x m degree square mentioned in equation 
98. If current practice were followed, the value of n i and hence N t , would 
depend on the following factors. 

(1) The error of representation E {Ag) n of an 0 x n° square, as defined in 
section D to the Appendix. This is a measure of the variability of A g 
within a square of a given size. Individual values of E (Ag> n are well 
known to depend on topographical variations in the case of the gravity 
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anomaly but the magnitude of E (Ag} n can in general be assumed to be 
representative of a given value of n for purposes of statistical estimation. 

(2) Stokes’ function f (0) should vary linearly over the region. If a (0, a ) 
system of coordinates were used, /i = sin 0 and it is convenient to 
define 

F(i /») = f (i/0 sin 0 (101) 

which is more stable than f (0) for small \p. 

(3) No correlation should exist between the variations in F (0) and Ag from 
the value adopted for the representation of the square. 

O O 

Consider in the first instance the representation of A g for the n x n square 
as afforded by the mean value Ag for the square, situated at the square center 
at which point the value of F(0) is F (0). If each n° x n° area were subdivided 
into N m° x m° equal area sub-divisions, let the individual values of the gravity 
anomaly and F(i p) be related to Ag and F (0) by the relations 

Ag t = Ag + c g . and F(\p k ) = F(0) + c^.. 

The total contribution to the final integral by the n° x n° area is given by 

N 

t = R m 2 + + <W 

i = l 

N N N 

= Kn 2 Ag F (0) +m 2 K F 0/T) ^ c gi +m 2 KAi 22 * ” 2k 21 C gi C \pi (1^^) 

i=l i=l i=l 
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where 

N = (n/m) 2 (103) 

The first term in equation 102 is the contribution due to the adoption of the 
area mean at the center of the n° x n° square while the second and third are 
zero by definition. The final term will tend to zero if there is no correlation 
between Ag and F (\p) as mentioned above. This possibility can be lessened further 
by restricting n so that variations in f (i/>) are linear over the area to the desired 
precision. In the present case it is desired to keep e N down to o {e 3 h d )= o {+5 cm). 
Hence the departures of Ffy) from linearity over the n° x n° should not exceed 
e 3 Ffy). The magnitude of variations in f(0) and Ffy) are functions of \p. Table 1 
gives the relationship between the square size n and \p which satisfy the linearity 
relation. 

The use of Stokes' function f(0) to evaluate equation 98 for all 0 would in- 
volve approximately two million summations if the above limits were adhered to 
and the effect of representation errors from Stokes' function were to be kept 
below the requisite magnitudes. The function F(0) defined by equation 101 is 
more stable for small \p but less economical to use than f(0) for large i// . N t 
can be reduced by a factor of 3 if F(0) were used to evaluate equation 98 when 
0 < 1° while i(xp) were used for all other 0 . This calls for the use of data de- 
fined on a local coordinate system (4> , a) for small 0 instead of the more 
universally applicable (0 c , system. Consequently, the definition of the former 
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data set from the latter must precede local computations which is why it is 
preferable to restrict such conversions to as small a region as possible. 

The effect of variations in Ag with position within the n° x n° square, on e N 
and the consequent representation errors, are best studied by analyzing the 
statistical characteristics of A g. The gravity anomaly as determined at any 
point on the Earth's surface is based on the following data. 

1. Observed gravity. 

2. Geocentric position of the gravity station. 

3. The difference in potential AW with reference to the "geoid." 

If Stokes' integral is to be solved by quadratures, it is relevant to investi- 
gate the errors which arise in the computed value of N f due to the representa- 
tion of a finite element of surface area by a single gravity anomaly. 

A useful function for the study of e N is the error of representation E {Ag) n 
for an n° x n° square (e.g., de Graaff Hunter 1935; Hirvonen 1956; Molodenskii 
et al 1962). More details of this important function and numerical magnitudes 
are given in section D of the Appendix. An empirical formula which describes 
the behaviour of this function is 


E(Ag) n = i Cj /n (104) 

A value of C : which fits most modern data in regions where the topographical 
gradients are small is 
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C x = 12 for 1/4° < n < 5° 

for n in degrees and E {Ag} n in mgal. For < 1/4°, a better definition of 
E{A g} n is obtained if the relation 

E{Ag} n = ±C 2 n (105) 

where C 2 ? 3X 10 under the same set of conditions as before. 

The first problem to be looked into is the effect of representation errors on 
e N if E {A g> n is assumed to have random error characteristics. In such a case, 
any other determinations of the gravity anomaly field in a specified n° x n° area 
which is represented for computation purposes by an adopted value Ag 0 , would 
deviate from Ag Q exhibiting characteristics implicit in the normal distribution. 

If square sizes in excess of n = 1° are excluded as these violate the prescribed 
linearity requirements of Stokes' function, as illustrated in table 1, it is interesting 
to verify whether variations in Ag over small squares are dominated by the 
gravity anomaly gradients BA g/3u a , u a forming an orthogonal and isometric 
angular coordinate system in the local horizon. 

Let the smallest sub-division of relevance be an m° x m° square whose 
error of representation satisfies the relation 

E{Ag) m = o{e 3 Ag} = ± o{50 /xgal} (106) 
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If equation 105 were true for very small n, m = 0.002°. The number of such 
regions in a larger n° x n° area would be N defined in equation 103. The gravity 
anomaly Ag and F (0) in the larger n° x n° square could be represented by Taylor's 
series of the form 


and 


U a 

Ag = Ag 0 + — 



F (0) = F(0 ) +— - 


< ~d> FQ/ Q 
dui 


(107) 


(108) 


where the origin of the u a coordinate system is at the center of the n° x n° 
square, the subscript o referring to values at this same point. On restricting 
the value of n to those square sizes where 


< o{e 3 F(i//)} for all j > 2 (i.e., n < 0.5°), 

Bijj 

a 


it is required to verify whether 


— < o{e 3 Ag} forall j^2 (109) 


when n < 0.5°. From equation A49 and 107 



60 



64 


If the m° x m° grid were symmetrically distributed about the n° x n° square, 
it is easy to show that 


and 



i = 1 


( 111 ) 


Hence equation 110 reduces to a relation of the form 


E {A g} n = ±Cn 

which agrees with the observed characteristics of E{ Ag > n as described by 
equation 105 for small n. Practical experience however indicates that sub- 
stantial deviations occur from this simplistic model especially in regions of 
rugged topography when the square size has to be reduced to an unacceptably 
small area to meet these specifications of linearity. It can be concluded that 
gravity an omaly variations are linear for square sizes under 0.3°, for purposes 
of statistical estimation, all departures having the characteristics of local noise. 
The contribution t of an n° x n° square to N f and hence h ^ as obtained by the 
evaluation of Stokes’ integral by quadratures, using N m° x m° squares on the 
basis of equations 98 and 100 is 

N 

t = K m 2 2 ^ A g. F (</».) 

i = 1 
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The use of equations 107 and 108 gives 


t - Km 2 


NA^F ( 0 O ) + F ( 0 O ) 


3 A 


d u 


IN 


A 3 F (0) * , \ ' 3 F (0) (Big , \ * 2 

A ■ cos A' / u + 1 L± _ cos A / u 2 . 

60 30 a Z_j ai B 0 \Bu a Z_i ai 

i = 1 \ “ i = 1 


- - A — s in A' y -1 u u \ + - F (0 ) - A g u 2 . 
Bu a U 2 7 2 3 u 2 4- 01 

i-l / ai=l 


o {e 3 t} 


= K n 2 


A g 0 F ^o) 


I F ( 0 Q ) [ — A B + £_Al\ + B F ( 0 ) B A g 


Bu 2 


B u 2 


B 0 3 u 


cos A' 


+ o {e 3 t} 


( 112 ) 


on using the results at 111, K being defined by equation 99, while A* has the same 
significance as at A 21. The magnitude of the contribution of the third term is 
governed by that of 3F(0)/30 which is two orders of magnitude smaller than the 
first for n > 0.1°, as can be seen from table 1. As the square meanAg n is given 
by 




1 3 j A g 
N Bui 


j ! 



(113) 


on using equation 107, it can be seen that the second term in equation 112 is also 
taken into consideration if Ag n were adopted instead of Ag o when representing 
an n° x n° square for the evaluation of Stokes' integral by quadratures. The terms 
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involving products of the type (3 F( 1 /' )/9uJ x ( 3Ag/ ?Ju a ) can be considered to have 
the characteristics of accidental errors only if no correlation existed on a large 
scale between variations in F('Z') and Ag. While such terms will make contribu- 
tions of significance when t in equation 112 is computed from a single n° x n° 
square represented by Ag^ instead of N m° x m° squares, the error is unlikely to 
have a regionality in excess of 1°. The use of typical values for the case when 
n = 0.1° indicates that the total contribution of this product term is less than 
0.1 mgal. The figures in column 2 of table 2 show that such errors will not 
affect the final results to ±1 cm even if the signs of 3F(0)/3<// and 3 Ag/3 u 

a 

were to hold the same sign over a 1° x 1° area. 

The above discussion may be summarized as follows: (a) The use of the 
mean value based on an evenly distributed sample gives a better representation 
than a single value when evaluating equation 112, the improvement being a func- 
tion of the number of the sample and the moment of the distribution of gravity 
stations about the square center, (b) The nature of the gravity field is such that 
any residual error due to the use of a tenth degree grid instead of smaller sub- 
divisions in non- mountainous regions has an effect less than 1 cm on the final 
value of N f . Also see section D of the Appendix. 

On adopting the basic square sizes specified in table 1, the remaining error 
characteristics can in the first instance be treated as random. The error e 

t a 

in t due to the error of representation E { Ag } m of the N constituent m° x m° 
areas is given by the addition law as 
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E {A g}J 


1/2 



(114) 


Over a limited n° x n° region F(</») can be represented by the first two terms of 
a Taylor's series when 


F (0) 2 - NF + 2 F (</>„) iLp 2Z d ^ + ° 
i^l i“ i 


{(HI 0 ) 



= N F (i// 0 ) 2 [l + o {10“ 3 }]- 


(115) 


as 



- 0 . 


On using equation 103 and noting that |F(<//) | < 2.5, 


where 


e ta = + o {K* m n E (A g) m (l+o {10 3 })} 


K' = 4 x 10“ 2 


(116) 


The accumulated accidental error e Na in N f is given by 


where 


Na 


= ± O 



± O (K" m E {A g }J 


(117) 


K" = 10. 
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Column 3 of table 2 gives estimates of e Na for various values of m in 
column 1 which represent the basic grid on which gravity stations are located to 
define the global gravity field for the evaluation of Stokes* integral by quadratures. 

It is also necessary to estimate the effect of an error e in evaluating t 
which retains its sign over an V x 1° area. For purposes of estimation, 
assume e ts to retain its numerical magnitude and sign over the larger area. 

The error e^ s in the larger block is obtained from equation 112 as 


where 



'ts 


= ±r 


e Ag being the systematic error in the gravity anomaly representing the n° x n 
square, K* having the same definition as in equation 116. The total contribution 
e Na to N f is estimated as in the case of equation 117 to be 


e Ns - ± ° (K"4 e Ag } (H8) 

Column 2 of table 2 gives estimates of e Ag for various values of l , specified in 
column 1, which ensure that e Ns = o{ e 3 h } . 

The following conclusions can be drawn concerning the evaluation of Stokes’ 
integral by quadratures. 
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1. The most critical factor is the departures from linearity of F(0), defined 
by equation 101. The use of F(< p) in practice is expensive as the (0, a) 
system of coordinates and not the (<P c , X) is used in computations, re- 
quiring the utilization of ring techniques which are less economic for 
computer use than the geographical square system. This is true even 
when use is restricted to those regions where <Jj < 3° and f(> p) is unstable. 
Computations with Stokes’ function in such regions calls for a finer 
sub-division in representing the inner zone gravity field on the lines 
described in table 1 as f(0 ) = o {10 3 } for \p = o {0.1° } while it is of 
order 10 2 for 4 1 - o {1°}. 

If this were not done, the terms ignored in equation 115 could be 
as large as the magnitude of those considered. Further, K’ in equation 
116 could in such a case be 2 - 3 orders of magnitude larger. Thus the 
four tenth degree squares within 0.1° of the point of computation would 
contribute ± o {0.3 cm} toward e N while the 100 tenth degree squares 
within 0.5 degrees would give rise to a further ± o {1 cm } due to 
departures from linearity (of order e 2 f(0)). 

2. In view of the limited errors introduced into the result, it can be con- 
cluded that a global gravity field based on definition at corners of a 
0.1° x 0.1° may be adequate for the evaluation of Stokes’ integral by 
quadratures to order e 3 h d (± 5 cm) if correct computing procedures 
were adopted and the gravity anomalies were free from systematic 

errors over large extents as specified in column 2 of table 2. 
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3. It is desirable that a procedure similar to Rice's circular ring method 
(Rice 1952) be used to compute the inner zones when F (0) should be 
used instead of f( { P) to circumvent the instability of the latter when 

\p < 1 °. This instability is not a consequence of variations in Ag and 
an adequate gravity field could be interpolated from the 0.1° x 0.1° 
grid without introducing significant error in the final result for studies 
of sea surface topography where effects of very short wave lengths are 
of no concern. 

4. The observation that interpolated values are not necessarily inferior to 
measured gravity anomalies was also made by Soviet geodesists 
(Molodenskii et al 1962, p. 172). The writer's own experience is that 
the extension of the gravity anomaly field represented by values on a 
uniform grid, to smaller subdivision in undisturbed regions, is stable 
without significant loss of accuracy (Mather 1967, p. 134). Thus if a 
0.05° x 0.05° grid were obtained by interpolation from a uniform 0.1° 

x 0.1° grid on which E{Ag} 0 = ±2.5 mgal, then E{Ag} Q Os o is held 

at this same value for the interpolated values, instead of the ±1.5 mgal 
estimated from equation 105. Thus the use of a 0.2° x 0.2° grid in lieu 
of the tenth degree grid as the fundamental basis of observations would 
result in e N = ± o { 6 cm }. 

5. While considerable laxity can be tolerated in the accuracy with which a 
reading represents a basic (i.e., tenth degree) square, the effect of 
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systematic error which retain the same sign over considerable extents 
has to be carefully watched. Table 2 shows that errors of 0.1 mgal 
which hold their magnitude and sign over 500 km can affect the com- 
puted value in excess of the specified limits of error. This type of 
error can be due to one of three causes. 

(a) Errors in the global gravity network controlling the gravity values 
used in the computations. 

(b) Loss of accuracy in unifying the elevation datums in relation to a 
globally acceptable "geoid." 

(c) Lack of precision in the geodetic coordinates used to compute the 
gravity anomaly as a consequence of regional datums not being re- 
lated to the geocenter as described in section 2.2. 

Thus IGSN 1971 can only be considered adequate in controlling the gravity 
fields in solutions to order e 3 if errors in defining individual station values in 
the net were uncorrelated and the interval between stations was not in excess of 
200 km. Neither of these criteria are likely to be met. On the other hand there 
are no limitations to present day meteorology which would inhibit the establish- 
ment of a global net with sufficient precision. Absolute station accuracies could 
be held at ±50 fj. gal as discussed earlier in this section using techniques similar 
to those used by Sakuma (1971). 

The unification of the elevation datums is equivalent to defining the geoid 
to a degree which has not been achieved as yet if the order of accuracy implicit 
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in this study is to be realized. The first stage in such a definition would be the 
adoption of a common epoch to which all so called "Mean Sea Level" datums are 
reduced. The principles underlying the establishment of such datums for a re- 
stricted region are well known and will not be discussed. The second stage calls 
for the definition of the sea surface topography and its departures with respect 
to a level surface which are stationary over long periods of time. The solution 
of the boundary value problem to order e 3 requires that errors of long wave- 
length in the definition of geopotential be kept to o {0.15 kgal m). This could be 
achieved without difficulty if the ocean surface were an equipotential, on allow- 
ing for tidal and meteorological variations. Unfortunately, the comparison of 
tide gauge readings with the results of geodetic levelling have indicated the 
existence of stationary departures of the sea surface from an equipotential as 
reported in section 1, both in the United States and Australia. As current prac- 
tice refers differences of geopotential to the sea surface instead of the equipo- 
tential, it becomes necessary to look into the effect such a procedure has on the 
computation of geoid heights with an accuracy of ±10 cm, which in turn calls for 
errors of less t han ±15 cm in the definition of the geoid as a datum for elevations. 

4,4 The Role of Satellite Solutions for the Gravity Field in Solutions of the 

Boundary Val ue Problem to Order e 3 

The characteristics of the Earth’s gravity field can be established by two 
different techniques. 

(a) The measurement of gravity at discrete points at or near the Earths surface. 
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(b) The determination of the disturbing potential from the analysis of the 
orbital perturbations of near Earth satellites. 

Solutions at (b) are interpreted in terms of spherical harmonic coefficients 
which can then be downward continued to the surface of the Earth with minimum 
mathematical complications. As Laplace's equation is satisfied at finite eleva- 
tions exterior to the Earth's surface, the disturbing potential V ' of the solid 

d 

Earth satisfies equation 41 which can be written as 


00 



n = 0 



n ji 1 


(41) 


where it is customary to express A in the form 

n 


A 

n 


n 



m = 0 


c 

nm 


(119) 


^nm being defined by equation 25. The absence of the first degree harmonic 
places the center of the reference ellipsoid at the center of mass of the solid 
Earth. 

The disturbing potential V ds as used in the analysis of orbital perturbations, 
is defined as that which causes the geopotential to deviate from that of a sphere 
with the same mass as the Earth. A symmetrical mass distribution is also im- 
plied when referring perturbations to the model adopted for central force motion. 
It is this potential whose derivatives define Lagrange's equation of planetary 
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motion (e.g., Kaula 1966, p. 29). As the gravitational effect of the atmosphere 
is estimated at less than 10 ^ gal at 30 km elevation, with the effect falling off 
rapidly with increase of elevation (IAG 1970, p. 72), it can be assumed that 

V d = V' + o{ e 3v d } 

at orbital elevations. The term of zero degree V d0 in V d ' has no effect on orbi- 
tal perturbations though its numerical magnitude could have a scaling effect on 
the orbital parameters used in the evaluation of the coefficients at 119. A further 
difference between V ds and V d is due to the ellipsoidal reference model used in 
defining the latter as described in section 2.2 in contrast to the spherical model 
used in obtaining V ds . If the effect on the gravitational potential is V de , then 


V 


ds 


= V' - V 

d V d0 


+ v. 


de 


( 120 ) 


On taking these factors into account when evaluating the coefficients C in 

a nm 

equation 119, the height anomaly h d at the surface of the Earth is given from 
equation 33, 41 and 119 as 
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h d >, n^l (121) 


where R, V a and y refer to values at the relevant point at the surface of the 
Earth. The infinite series must by definition converge to the limit specified by 
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equation 121. The evaluation of the coefficients C a nm defined in equation 25, by 
the analysis of orbital perturbations, is dominated by two effects. 

(1) The damping effect of the term (a/R) n as R > a (e.g., see Mather 1971c, 
p. 67). Consequently the coefficients of degree n less than some limit- 
ing value n s are well determined, n s being a function of the orbital 
elevation of the satellite. 

(2) The effect of resonance between the values of the set {n, m} and the 
orbital period. This causes certain coefficients which by themselves, 
make no contribution of significance towards the representation of the 
Earth's gravitational field, to have marked effects on the perturbations 
of those orbits with sympathetic parameters. As a consequence, all 
orbits are sensitive to certain resonant harmonics whose identity can 
be predicted from the orbital elements (e.g. Wagner 1967). 

At first glance it would appear that a very large number of satellites in a 
variety of orbits would afford a means for the complete determination of the 
Earth's gravitational field. The costs involved make such evaluation unlikely due 
to masking effects which make it difficult to separate some of the resonant terms 
unless adequate variations were available in the orbital inclination. Serious 
thought should be given to the role concepts of resonance should play in solutions 
of high resolution for the Earth's gravitational field at the surface of the Earth 
as it is most likely that only a limited number of satellites will be available for 
the task. These accurate determinations will suffer from a loss of resolution if 
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not restricted to a limited interval of time. If such is the case, it may be 
preferable to treat higher degree resonant effects as sources of orbital pertur- 
bations rather than signals from the gravity field which could be meaningfully 
translated into representations at .the surface of the Earth. 

The gravity anomaly at the surface of the Earth is obtained from equations 

41, 45 and 119 as 


A kM 

Ag = > (n - 

R 2 L-x 

n = 0 
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+ i g C 2 + o { e 3 Ag} , (122) 


n i 1 

The comparison of the values of Ag computed from the C a nm determined from 
the analysis of orbital perturbations with those established from surface gravity 
measurement on allowing for equation 120, provide an index of the success with 
which a truncated spherical harmonic series (i.e., n< n.) can represent the 
gravity anomaly at the surface of the Earth. This could be extended further by 
incorporating those harmonics in the range n s < n < n^ from surface gravity to 


enhance the representation provided by the spherical harmonic series, thereby 


increasing the range of the power spectrum and reducing the residuals on com- 
parison of deduced and observed values of the gravity anomaly at the surface of 
the Earth. Such concepts assume that those C a nm for n < n s as determined from 
orbital perturbations were free from error as were the values of Ag at the sur- 
face of the Earth. It also has the advantage that errors in the framework con- 
trolling the gravity datum at the Earth's surface, which is established with 


73 



77 


difficulty, can be eliminated or at least minimized in the representation of the 
global field. 

This procedure is not always strictly followed in common practice when 
general adjustment techniques are used to minimize residuals without holding 
any of the quantities fixed. 

The variance of gravity anomalies at the surface of the Earth is approxi- 
mately 1200 mgal . Solutions to (20, 20) absorb over 90% of the power inherent 
in the representation (Lerch et al 1972, p. A12). Thus 

M {(A g Q _ A g s ) 2 } = 100 mgal 2 , (123) 

where the subscripts o and s refer to terrestrial and satellite determined values 
respectively. The absorption of this balance 10% of the power spectrum is 
likely to require a great increase in the number of terms though some of this 
residual is due to deficiencies in the surface gravity data. From a study of the 
error of representation, given in section D of the appendix, a (20, 20) solution 
can be considered to be equivalent to a representation on an 0.5 degree grid only 
if the comparisons represented at 123 were with individual gravity values. This 
is not the case, as the surface gravity values were in the form of five degree 
area means. The conclusion which can be drawn is that the (20, 20) representa- 
tion is equivalent to a global 5° x 5° coverage with 5-6 readings per square and 
zero moment of distribution about the square center. 
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The use of spherical harmonic representations of the gravity field to achieve 
the definition of the gravity field required in conclusion 2 to section 4.3 (i.e., a 
tenth degree grid) would require that the former absorbs all but 9 mgal 2 of the 
power spectrum on comparison with individual values. This will be equivalent 
to absorbing all but 0.1 mgal 2 of the power spectrum on comparison with ade- 
quately co mp uted one degree area means, each based on 100 evenly spaced values 
with zero moment of distribution about the square center. The latter would in- 
volve an al ysis up to degree 180 (over 3 x 10 4 coefficients). It has yet to be 
established whether such refined determinations of the gravity field are possible 
by satellite to satellite tracking of low altitude satellites. 

The requirements for a complete solution of the geodetic boundary value 
problem to order e 3 (i.e., ±5 cm in h d ) is a gravity field representation based 
on at least a tenth degree grid. This is equivalent to a spherical 
harmonic representation where n = 1800, involving o { 3 x 10 6 }"terms, which is 
not si gnif icantly different from N t in table 1 . The use of such functions can 
therefore be justified in this case only if the amount of surface gravity informa- 
tion on the tenth degree grid were significantly low. It is unlikely that any 
favorable claim can be made at the present time regarding the achievement of 
this degree of resolution from the study of orbital perturbations. It would there- 
fore appear that satellite determinations of the gravity field could well be inferior 
for the complete solution of the geodetic boundary value problem to ±5 cm if 
(a) surface gravity data were available globally on a tenth degree grid; and 
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(b) systematic errors in the gravity anomalies were held to below ±50^ gal. 
The low degree harmonics from orbital perturbations could however play a 
significant role in such solutions when (a) is true but not (b). The three major 
sources of systematic error in gravity anomalies which have long wavelength 
are given in note 5 to section 4.3. While (c) is likely to be resolved with minimal 
difficulty, systematic errors at (b) are complex primarily as a consequence of 
possible stationary departures of the sea surface from an equipotential. If the 
gravity anomalies have been corrected for effects at (b) and (c), any residual 
long wave discrepancies between surface gravity data based on adequate samples 
and the low degree harmonics obtained from the analysis of orbital perturbations 
and with the required precision, should provide an effective check on the 
systematic error propagation of the type at (a) in the note mentioned above. 

The results obtained from the analysis of the orbital perturbations of satel- 
lites in near Earth orbits are unlikely therefore to provide the representation of 
the gravity field which is required for a complete solution of the geodetic boundary 
value problem to order e 3 . The determination of the low degree harmonics in 
this representation with adequate precision will however be invaluable in resolv- 
ing the systematic errors in the global gravity standardization network described 
in note 5 to section 4.3. 

4.5 Departures of the Sea Surface Topography from an Equipotential Surface 
Until recently, no attempt has been made to study the departure of the sur- 
face of the oceans from a level surface. The existence of such departures has 
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been established on comparing the results of geodetic levelling with tide gauge 
readings. These departures which will be called stationary, in order that they 
could be differentiated from short term effects due to winds, other meteorological 
factors as well as the short period distortions on the geops due to tidal effects. 
The use of satellite altimeters provides a means of determining the instantaneous 
geocentric position of those features of the sea surface with wavelengths in ex- 
cess of ^ km. -1 = 200 for the proposed GEOS-C mission (Weiffenbach 1972, 
p. 1-1). The stationary departures can be obtained by allowing for the effect of 
tides and meteorological conditions, on differencing equivalent position vectors 
to the sea surface and the geoid. 

As only features with wavelengths in excess of F km are being studied, it 
is possible to use a truncated version of equation 121 to obtain the required 
definition of the geoid even to order e 3 h d . Over oceanic regions, the telluroid 
coincides with the geoid and the elevation N of the latter above the ellipsoid is 
given by 


N = hj, if AW = 0 ( 124 ) 

If the gravity field were represented by a global set of gravity anomalies, N 
could be obtained from the set of equations summarized in Section 4.1. Alter- 
nately, these anomalies could be analyzed for the equivalent harmonic coefficients 
using equation 122 and the values of N in oceanic areas obtained from equation 
121. From the discussion in section 4.4, the representation should absorb all 
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but 0.1 mgal of the power spectrum on comparison with one degree square means 
compiled from 100 values spaced on a tenth degree grid, if wavelengths in excess 
of <£ km in N is to be defined to e 3 N. This should give an accuracy of ±10 cm in 
N on the basis of the results in table 2, which is a desirable goal in the definition 
of the sea topography (Williamstown Report 1969, 3-2). 

Consider the use of equation 121. The harmonic coefficients could be ob- 
tained from surface gravity on controlling gravity standardization network 
errors with low degree harmonics determined from orbital perturbation analysis 
of adequate precision. In practice it is likely that the distribution of surface 
gravity information will continue to be non-uniform. It is therefore relevant to 
designate a desirable form in which the gravity data should be used in the 
analysis for harmonic coefficients. A global representation on a tenth degree 
(10 km) grid has an error of representation of approximately ±3 mgal, resulting 
in an accuracy of ±10 cm in N if the data is free from large scale systematic 
errors. A study of equation 116 indicates that if m = 0.1°, the precision required 
in the mean value for a n° x n° area to maintain the specified accuracy e N in N 
is not E{Ag) n but (E {Ag} m x m/n), all other things being equal. Thus the 
equivalent precision required from a 1° x 1° square mean is approximately 
±0.3 mgal'. Such a mean can be computed only if 

(1) 100 values spaced on a tenth degree grid are used in its evaluation; and 

(2) the moment of distribution of the gravity stations about the square 
center tends to zero. 
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This does not mean that each one degree square should contain 100 gravity 
readings on a tenth degree grid. It is well known that gravity anomaly values 
can be predicted under carefully controlled conditions such that the prediction 
error did not exceed the error of representation (e.g., Molodenskii et al 1962, 
p. 172; Mather 1967, p. 134). The exact technique to be used for this purpose is 
a matter for debate. In practice, the writer has found that practical and not 
theoretical considerations predominate in the choice of a particular method. 

Any commonly used interpolation routine will give the desired accuracy 
provided 

(a) sufficient data were available to avoid predictions based on readings 
which were not in the immediate vicinity of the point; and 

(b) the correlation of gravity anomalies with elevation over limited regions 
were allowed for. 

For example, an evenly space 50% coverage of a 1° x 1° square (i.e., 50 
readings) should give the required accuracy in the area mean if the latter were 
computed from 100 evenly spaced values with zero moment of distribution about 
the square center and the above requirements were met. Tests carried out for 
non- mountainou s regions in Australia with considerable gravity variation, indi- 
cated that a 20% representation, again evenly spaced, could provide interpolated 
values whose error would be double that for E { } Q 1 (ibid, p. 133). In such a 
case, E{Ag } j = ± o {0.5 mgal). This figure falls off markedly if the moment of 
distribution about the square center did not approach zero. 
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The following conclusions can be drawn about the preparation of area means 
prior to harmonic analysis in regions which are incompletely represented by 
surface gravity data. 

1. Values should be predicted from available observations represented on 
a tenth degree (10 km) grid using any reasonable interpolation routine 
or collocation techniques, and allowing for height correlation as well as 
the deviation of gravity station elevation from the mean elevation of the 
region it is intended to represent. 

2. If a network of gravity stations were being planned, the stations should 
be cited such that the distance over which interpolations are made 
should be as small as practicable to avoid systematic effects. 

3. The quality of the area mean is more dependent on the nature of dis- 
tribution of the gravity stations about the square center, rather than 

the number of readings available. This is characterized by the moments 
M a of the gravity station distribution defined by 

N 

M a = ^ d U ai ( 125 ) 

i = 1 

where du ai are the coordinate displacements of the i-th gravity station 
from the square center. More research is necessary into the role M a 
should play in setting up observation equations for the determination of 
the harmonic coefficients. 
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The error e^ in A g as computed from equation 122 due to an error e c 

B n n 

in C n is given by 


e A» = ± o{y(n - 1) e c n }, 

s n 

where e r is the r.m.s. of the sum of the variances of the 2n + 1 coefficients of 

'~n 

degree n. The analysis of harmonic solutions of this type indicate that the mag- 
nitude of the average variance of coefficients of degree n are essentially constant 
(say cr 2 ) for solutions up to degree 12 (Lerch et al. 1972, p. 21) while de- 
partures of individual variances from this mean, fall within a maximum amplitude 
cT_ ax . On assuming sinusoidal characteristics for the deviation e. of individual 
standard standard deviations from o - , the total variance of terms of degree n is 


2n + l 


2 n = o (cr + e.) 2 = o«j(2n + 1) cr 2 + ^ 

= „ {(2n + 1) (<r 2 + Kcr^)} = o {(2n + 1) ej)| 


(126) 


where e c is a constant whose magnitude is approximately 2 x 10 -8 for solutions 
obtained at the present time. The extension of these observations seem to indi- 
cate that e A is a function of n 2 . On the other hand, the results in column 2 of 
table 2 indicate that larger errors could be tolerated in the higher degree har- 
monics without significantly worsening the results if equation 90 were used in 
the computation. The required accuracy for those of low degree is about 5 parts 
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in 10 4 if each is treated as an isolated error source. It is difficult to estimate 
composite effects in solutions to very high degree in the absence of solution 
characteristics. It could be assumed that an adequate algorithm will result in 
the harmonic representation having error accumulation patterns similar to those 
of the original data, provided the latter were free from systematic error. 

If surface gravity data were used to determine the geoidal slopes with 
wavelengths longer than 200 km, it would therefore be necessary to compute 
1° x 1° (100 km) area means from 100 evenly spaced values on a 0.1° x 0.1° 

(10 km) grid in non- mountainous areas, such that the error of representation of 
the area mean is ±0.3 mgal. This would ensure ±10 cm accuracy in the com- 
puted result. The analysis of such a data set for the appropriate coefficients 
using equation 122, followed by the evaluation of h d from equation 121, should 
give the required result. The existence of such a data set could also be used to 
give the same result through equations at section 4.1. In both cases it is ex- 
tremely desirable to verify the correctness of the low degree harmonics against 
satellite determined values of adequate precision, to ensure that the results are 
free from systematic errors in the compilation of the global elevation and gravity 
datums. 

Notes : 

(1) It should be pointed out that it is quite valid to use the truncated spherical 
harmonic series in equation 121 for the evaluation of the characteristics 
of the geoid with wavelengths in excess of a certain minimum value, 
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provided such values in themselves are capable of meaningful interpre- 
tation. As this information is to be used in conjunction with altimeter 
data which can only evaluate similar characteristics of the sea surface, 
it is relevant to attempt the definitions of the long wave characteristics 
of geoid to ±10 cm, noting that such evaluations could deviate from the 
true stationary geoid over oceans by up to ±5 m. 

(2) The development given above has only dealt with the techniques for the 
solution of the sea surface topography using determinations of the gravity 
field at the surface of the Earth. Satellite techniques which have been 
proposed for reaching similar goals (Williamstown Report 1969, 2-20 - 
2-24) have not been considered as they fall beyond the scope of the 
present development. The equations in section 4.4 are of relevance 
however when formulating a solution to the problem in this case as 
well. 

4.6 A Note on the Determination of Gravity Anomalies at the Surface of the Earth 
The establishment of the gravity field at the surface of the Earth for de- 
terminations of sea surface topography with a resolution at the ±5 cm level does 
not require that individual gravity determinations are consequences of techniques 
achieving accuracies of better than ±50 jugal and equivalent station elevations to 
±15 cm at each point. Instead what is required is the control of the propagation 
of systematic error due to those sources with long wavelength, to values below 
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these exacting limits, on the basis of equation 118, as illustrated in column 2 of 
table 2. 

The factors which influence such errors are the establishment of latitude on 
a geocentric datum, elevation and gravity such that these systematic error limits 
are not exceeded. The error e , in latitude is discussed in section 2.2. It is 
current practice to compute normal gravity from the value of 4> referred to the 
local geodetic datum. If e^ =2 arcsec at 0 = 45°, the resulting systematic error 
in Ag = ± 0 {50 [i gal}. It follows that the application of orientation vector cor- 
rections from any of the more recent satellite solutions for geocentric position 
prior to the computation of normal gravity, will ensure that this source of 
systematic error is eliminated. 

The effect of elevation errors e h in the gravity station elevations used in the 
computation of the gravity anomaly is not straightforward. Errors approaching 
±50 Aigal are obtained in Ag if e h = ± o (15 cm). Such a specification is at the 
noise level of internal errors in large first order regional geodetic level net- 
works. As pointed out in section 4.3, the essential requirement is the control of 
systematic errors with long wavelength when establishing the global datum for 
elevations. This would call for a consideration of 

(a) the time dependent variations in "Mean Sea Level"; and 

(b) the stationary departures of the sea surface from the equipotential 
surface adopted as the geoid. 
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Techniques for the estimation of the former constitute areas of regular re- 
search in oceanography. The effects at (b) need evaluation only at those points 
on the coast which have been used to define the sea level datum for elevations. 
The geocentric position of a tide gauge on each of these reference datums could 
be established in the future by means of a connection to a suitable laser ranging 
station which forms part of a global network. The elevation of the sea surface 
can only be determined if the geoid height at this point were known. There is 
little choice but to iterate between improvements in the elevation datum and the 
determination of the geoid to obtain a solution of adequate accuracy, a procedure 
which could be quite expensive as there may be difficulties in making the solu- 
tions converge, as illustrated below. 

Elevation errors of considerably larger magnitude can be tolerated in sta- 
tion elevations provided they are purely local in character. It should be noted 
that an error of 1 m in the elevation is approximately equivalent to 0.3 mgal in 
the gravity anomaly, which in turn can have an error of representation of ±3 
mgal in the context of the global grid discussed earlier. 

The use of accurately determined low degree harmonics of the Earth's 
gravitational field from the analysis of orbital perturbations for the verification 
of the global gravity standardization net will be successful only if the errors in 
the establishment of the global elevation datum have been satisfactorily resolved. 
The sea level datums in current use cannot be considered to be compatible for 
the purpose of solving the geodetic boundary value problem to order e 3 , as no 
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serious attempt has been made to define the term "geoid" on a global basis to 
±15 cm. Elevations with respect to local determinations of the sea surface can 
be considered to be elevations above geoid only if 

(1) corrections were applied for the epoch of determination; and 

(2) the stationary departures of the sea surface from the equipotential were 
allowed for. 

The latter is difficult to accomplish in advance of a geoid determination to ±15 
cm unless all the tide gauges are linked by a single network of geodetic levelling. 
While such connections would be feasible for the American continents as a unit 
or Africa, Asia and Europe as a second entity, a global connection cannot be ef- 
fected to achieve this end. If errors on this count averaged at ±1 m, causing 
effects in the gravity anomaly of approximately ±0.3 mgal with wavelengths of 
1000 km, the accuracy of the computed value of N is estimated at ±15 cm. In 
such a case, the error in the determination of stationary departures of the sea 
surface from an equipotential can also be determined to ±15 cm, assu ming that 
the geocentric positions of the tide gauges defining the elevation datum are es- 
tablished with this same order of accuracy either from laser ranging tec hni ques 
or from satellite altimetry. The systematic errors in the gravity anomalies due 
to the revised height datum are almost an order of magnitude smaller and hence 
fall within the required limits of precision for a solution of the boundary value 
problem to order e 3 . 
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If, on the other hand, the dominant stationary characteristics of the sea sur- 
face topography had twice the wavelength and magnitude as in the above case, the 
error in the computed value of N is estimated at ±60 cm in the initial iteration. 
Consequently two iterations are necessary to ensure the definition of the geoid 
to order e 3 h . 

. d 

The use of collocation techniques in defining the unsurveyed portions of the 
Earth’s gravity field is outside the scope of the present development. The ac- 
curacy of any predicted values are most likely to meet the criteria given in note 
4 to section 4.3 if based on a minimum of four equidistant values, each pair of 
which subtends nearly equal angles at the point of prediction, and in regions 
where topographic variations are smooth. As pointed out in section 4.3, a net- 
work pre-planned in such a manner could be used to increase the gravity anomaly 
representation by a factor of 4: 1 in undisturbed regions without introducing 
significant error provided the gravity values at those points used to control the 
prediction, are substantially free from the sources of systematic error described 
above. Other criteria of significance are the following. 

(i) Predictions should be restricted to regions where the behavior of the 
gravity field is sufficiently well known so that the error of prediction is 
no greater than E {Ag}, as discussed in section 4.3. 

(ii) The prediction interval is small enough to permit the adequate repre- 
sentation of the correlation of gravity anomalies with elevation. 
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5. CONCLUSIONS 

The above development has defined formulae for the height anomalies h d at 
the surface of the Earth to order e 3 h d (±5 cm). h d is equivalent to the elevation 
N of the geoid above the ellipsoid of reference in ocean regions. The solution 
obtained for the boundary value problem in geodesy includes Stokes’ integral in 
circumstances where spherical harmonic expansions are resorted to only when 
physical validity exists for their use. The boundary value condition has been 
built around the potential of the solid Earth and oceans, excluding the atmosphere, 
to ensure the mathematical validity for the expression of the solution in terms of 
spherical harmonic expansions. Such a representation is desirable as it permits 
the ready incorporation of information regarding the Earth’s gravitational field, 
as obtained from the analysis of the orbital perturbations of near Earth satellites, 
in representations at the surface of the Earth. 

The solution referred to above exists only if the stratification of the 
atmosphere, assumed invariant with respect to the epoch in which the gravity 
field is determined, were known. A model has therefore to be defined for the 
Earth's atmosphere prior to effecting the solution, which is referred to the center 
of mass of the solid Earth. The latter can be related without difficulty to the 
geocenter if the mass distribution of the atmospheric model were known. This 
could be defined as a series of ellipsoidal shells (e.g., IAG 1970, p. 62). More 
correctly, the atmosphere can be considered to consist of such shells at alti- 
tudes greater than the maximum topographical elevation. At lower altitudes, 
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the shells are not complete, the excluded portions being related to the topography 
of the Earth. 

Contributions to the height anomaly h rf , as computed using the formulae 
listed in section 4.1, have three distinct orders of magnitude to o { e 3 h d ) . The 
major term is obtained by the use of gravity anomalies at the surface of the Earth 
in Stokes' integral which gives over 90% of the total magnitude. The second is 
due to departures of the topography in the local area from a plane, contributing 
less than 10% of the total magnitude in regions of rugged terrain. The third set 
of terms is of order e 2 h rf (±30 cm) and arise as a consequence of the ellipticity 
of the Earth, topographical gradients at the surface of the Earth and the conse- 
quences of the Stokesian assumptions. 

No direct solution is possible, as pointed out in section C of the Appendix 
and an iterative procedure, described in section 4.2, has to be resorted to. The 
representation of the gravity field which would ensure adequate accuracy in the 
evaluation of Stokes' integral by quadratures can be estimated from those char- 
acteristics of gravity anomalies embodied in the error of representation E{Ag}. 

It is estimated that the definition of this field at the surface of the Earth by values 
on a tenth degree (10 km) grid in non- mountainous regions (97% of the globe) would 
ensure that the accuracy of the value computed for the Stokesian terms was of 
order e 3 h d , provided no systematic errors with long wavelength existed in the 
data. Regions of rugged topography and disturbed areas close to the point of 
computation should be represented by square sizes whose error of representation 
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remains at ±3 mgal, if the accuracy of the final result is not to deteriorate by 
as much as a factor of 2. 

It is therefore essential that the following criteria are met in defining the 
global gravity anomaly field. 

(1) Stations comprising the global gravity standardization network should 
be established with absolute accuracies of ±50 /xgal, and at intervals 
which are not much in excess of 1000 km 

(2) The datum to which measurements of geopotential are referred (i.e., 
the "geoid") should be defined with an accuracy of ±15 cm. All "Mean 
Sea Level" datums should be reduced to the common epoch of the gravity 
measurements and the stationary departures of the sea surface from 
the equipotential corrected for prior to the computation of gravity 
anomalies . 

(3) All values of normal gravity should be computed using the equivalent 
latitude on a geocentered ellipsoid, correct to ±2 arcsee, rather than 
values on a regional geodetic datum. 

(4) Individual gravity stations should reflect the mean elevation of the 
region. It should be noted that comparatively large errors can be 
tolerated in individual gravity anomalies on the tenth degree grid 
provided they are random in character (i.e., purely local). 

Gravity observations need not be made at every point on the tenth degree 
(10 km) grid in non- mountainous areas. It is common experience that 


90 



94 


interpolation techniques can give predicted values without increasing the error 
of representation under carefully controlled conditions. This factor should be 
taken into consideration when planning any large scale sampling of the global 
gravity field. 

It is unlikely that numerical solutions aiming for an accuracy of ±10 cm in 
h d can be achieved without significant roles being played by techniques from 
satellite geodesy. The first is in the establishment of geocentric position at 
tide gauges monitoring the sea surface and hence the datum for geopotential 
differences, and therefore gravity station elevations. This information would 
provide the link between the results of geodetic levelling and the geoid on a global 
basis, on using the iterative procedure described in section 4.5 if necessary. 

The method so described could prove ineffectual if the stationary departures of 
the sea surface from the geoid are characterized by very long wavelengths and 
amplitudes in excess of 2m. 

The second role that could be played by techniques in satellite geodesy is 
the determination of the low degree harmonics of the gravity field from the 
analysis of orbital perturbations for the control of systematic errors with long 
wavelength in the gravity anomalies at the surface of the Earth. If the sources 
of error at (2) and (3) above are allowed for, such coefficients of adequate pre- 
cision could control systematic effects in observed gravity due to standardization 
network errors, especially in ocean areas, where stations in such a network 
may be widely spaced. Consequently, the absolute error in the gravity control 
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station, which should be less than ±50 /xgal, could be repeated in all gravity 
anomalies over a very large area. 

The determination of features of the geoid with wavelength in excess of 
200 km calls for the evaluation of equation 121 using a truncated harmonic series 
including all terms up to n = 180 (i.e., over 3 x 10 4 coefficients). Such de- 
terminations are meaningful in the context of quantifying the stationary de- 
partures of the sea surface from the geoid, as only wavelengths in excess of the 
figure given above can be resolved by altimetric techniques. The minimum surface 
gravity anomaly field necessary to obtain the resolution of such features to ±10 cm 
is estimated to be a 1° x 1° (100 km) grid which is represented by the area mean 
value computed from at least 100 equally spaced values with zero moment of 
distribution about the square center. The harmonic representation should be 
capable of absorbing all but 0.1 mgal 2 of the power spectrum of the gravity 
anomalies at the surface of the Earth as represented by correctly computed 
area mean values on the one degree grid. 

Satellite geodesy, as distinct from satellite altimetry, still has an important 
role to play in defining stationary departures of the sea surface from the geoid, 
even if satellite-to- satellite tracking and drag-free satellites do not play a major 
role in the definition of the global gravity field. A world-wide system of tracking 
stations for laser ranging to satellites could provide both the resolution of the 
datum for the measurement of geopotential, as well as an accurate determination 
of the low degree harmonics of the Earth’s gravity field. 
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There appear to be no long range obstacles that are likely to inhibit the 
definition of those characteristics of the geoid in ocean areas with wave lengths 
in excess of 200 km with a ±10 cm resolution. An evenly spaced sampling of 
the global gravity field on the line? described in section 4.3, and based on 
levelling and gravity control networks with suitably small systematic error 
c har acteristics, remains a necessary pre-requisite for a successful determina- 
tion from surface gravity data. Information controlled in this manner and located 
on a grid where the station spacing could be as large as 20 km in non- mountainous 
and undisturbed regions, and when used in conjunction with prediction methods 
which took elevation correlations into account, is likely to provide the desired 
resolution. 
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Preceding page blank 

~~~ APPENDIX 

A. Relations on the Ellipsoid 

The relation between geocentric and geodetic coordinates are readily avail- 
able in many texts on geometrical geodesy. These expressions do not, as a rule, 
consider the effects of the topography. Given an ellipsoid of revolution of equa- 
torial radius a and flattening f, the geocentric latitude 0 CO of a point P 0 on the 
reference ellipsoid, is related to its geodetic latitude by the formula 


S0 = <t> - <£ co 


tan 0 - tan 0 co 
1 + tan 0 tan 0 co 


As tan 4 > cq = (1 -' f) 2 tan 0, 8 0 = tan 80 + o{f 3 }. 

Thus 


80 = Z 1 ~ f _ 1 — L sin 0 c cos 0 c 

1 + 2 f sin 2 0 

C 

= f sin 20 c ' |l +-| f - 2 f sin 2 0 c + o{f 2 }J ( A_1 ) 

Further, if the elevation of the point P at the surface of the Earth, above P 0 
is h, as illustrated in figure A-l and if PGP 0 = §0 c , where G is the geocenter, 
assumed coincident with the center of the ellipsoid, it follows that 

sin S 0 c _ sin §0 

h R 

PRECEDING PAGE BLANK NOT FILMED 

A-l 
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where R is the spherical coordinate of P on a geocentric spherical system of 
coordinates (R, 4 > c , k). On using equation A-l, 

8 4> c = — f sin 2<£ c + o{f 3 } (A-2) 

R , 

Also, 


R = R q cos S0 c + h cos (8 4> - S0 c ) 


= R 0 + h + o { f 3 R } 


(A-3) 


where R Q is the distance to P from the geocenter G. For most practical 
purposes, 


R 0 = a [1 - f s in 2 0 c + o{f 2 }] , 

and the mean radius R m of the Earth ellipsoid is given by 

R m = 3 f + °^ f2 > • 

The combination of the above three relations gives 


R n = R 

0 m 


l+f ( ^ - s in 2 



and 


R = R (1 

m v 


+ c r) 


(A-4) 


(A-5) 


A- 2 
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where 



m 



+ o{f 2 } 


(A-6) 


The distance r between the surface element dS at Q(R, 0 c , A. ) and the point 
of computation at P(R p , <£ cp , A p ) can be related to the angle subtended by the 
geocentric radii GP (= R p = R + AR) and GQ (= R) at G, as illustrated in figure 
A-2. As GP and GQ lie in the plane of the meridians through P and Q respectively, 
the angle dA. between the meridian planes is given by 


d\ = X _ K 

p 


(A- 7) 


Thus 


0 = cos -1 [sin <b sin <t> + cos <±> cos cb cos dX] 

cp c cp ^ c 


(A- 8) 


without approximation. 


The Term x 3 /r 2 

The term x /r 2 in equation 38 is obtained from figure A-3 as 


x 3 R cos S - R p cos (\p + S) 

_ 


(A-9) 


where 


S .= (S0 - §0 c ) cos a a 


= f sin 2 4 > c cosa^ +o{f 2 } 

A-3 


(A-10) - 
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a being the azimuth of P from the element of surface dS at Q. 
is given by 

r 2 = R 2 + Rp - 2R R p cos \p = (R p - R) 2 + 4R R p s in 2 - 
= (AR) 2 +4R 2 sin 2 I«A(1 +c R + c Rp ), 

where c D is defined by equation A- 6 and AR is given by 

R 

AR = R m (Sp - C R> = R m f < Sin2 ~ Sin2 ^cp) + h p ' h + 


On defining r by the relation 


r n = 2R sin — 

0 m o r 


the expression for r 2 can be written as 


.2 _ _2 


= r 0 C 1 + C r> 


where 


AR' 

C r = ' 4 


o{f 2 } 


The expression for x 3 /r 2 should be put in the form 


X 3 1 

— =— (1 + 4 >) 
r 2 2R 


The distance r 


4> 


o {f 2 R} (A-ll) 


(A-12) 


(A-13) 


(A-14) 


(A-15) 


A- 4 
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to facilitate the recovery of Stokes' integral from equation 38. The structure of 
for purposes of numerical computations is dependent on the magnitude of A R 
and i p. The third term in equation A14 could exceed unity when h p >> h for limit- 
ing values of 0 . To avoid loss of generality, it is preferable to retain closed 
expressions at this stage of the development. Equation 15 could be rewritten as 

(A-16) 


2R x. 


4> 


- 1 


The use of equations A-6, A-9 and A-13 in equation A-16 gives 


4> 


a + c r) 

r 0 + C r) 


[(1 + c R ) cos S - (1 + c Rp ) cos (0 + 8)] - 1 


2R 2 


(1 + c r ) 


-1 


2 s in 2 0 + 2c R + 8 s in 0 - c Rp cos 0 + o{f 2 } 


-1 (A-17) 


as S = o{f}. 

Simplified working expressions for 4> are obtained by fixing maximum mag- 
nitudes for (AR/r 0 ) 2 . For example, if [(h p = h)/r 0 l 2 = o{f }, for small r 0 , 
c = o{f }. Hence, equation A-17 can be written as 


2R 


2 r 


4> 


2c d + 8 sin 0 - c„ cos 0 - 2c sin 2 — xp - o{f 2 } 

K Kp r 2 


(A-18a) 


Alternately, if c r = o {10 -1 } 


2R 2 


4> = 


(1 - c r + c 2 + o{f}) (2 c r + 8 sin 0 - c Rp cos 0 + o{f 2 }) 


2 sin 2 i 0 ^ (-1) 1 c * + ° 


{f 2 } 


(A-18b) 


i = l 


A- 5 
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It would be advisable to use equation A-17 for the evaluation of $ when c r 
takes larger values. In conclusion, $ is usually a small quantity, except under 
limiting conditions when i p and hence r is small, together with the elevation of 
P differing considerably from that of the surrounding topography. If terms of 
the order of the flattening were neglected, the contribution of this term arises 
only from elements with small <// and great differences in (h - h p ), when equation 
A-17 can be written as (Mather 1971b, p. 81) 

g 

$ = _!L [h - h -R (c -c 2 + o{f}) + o{fh}] (A-19) 

« p m v r r ' 

r o 

The expressions for $ given at equations A-17 to A-19 can be programmed 
without problems to the required order of accuracy. 

The Terms x a /r 2 

Let the angle between the line QP and the XjX 2 plane be x, as illustrated 
in figure A-3. Then PQG = 1/2 77- (x - S) and QPG = 1/277 - (<// - x + 8). The 
application of sine formula to triangle PQG gives 

_ R (1 + c„) R (1 + Cp ) 

r _ m v R' _ m ' R P 

sirup cos (i p - y + $) cos(x-S) 

Thus 

cos (0 - X + s ) (! + c Rp ) = cos (x - s ) (1 + C R ) (A-20) 

It follows from the use of equation A-20 on referral to figure A-3, that 


A-6 
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x r cos Y cos A' p c • • , 
a _ a- _ K s in 0 cos X , 

T — -r-. COS A , 

r 2 r 2 r 2 cos (’/' - X - S) 


where 


Aj' = o. a and A' = 1 ti - 


(A-21) 


being the azimuth of P from dS at Q. Equation A-20 also indicates that 


cos(X-S) , AR rr9S 
— — - 1 +^r— + oiP) 


cos (0 - x + S) R 


(A-22) 


On ignoring terms of order f in equation A-22, 


cos x = cos (0 - x) + o{f} 

or 

x = i 0 + 0 {f}, 

provided 0 is not of the same order of magnitude as S . Let 

cos X 

cos(0~y + S) + C * (A-23) 

where C x is usually a small quantity of order f except when 0 = o {§} . Also 
define x by the equation 


X = 



A- 7 


(A-24) 
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The use of equations A-22 and A-23 gives 

cos ^-^-0 - 8 — < 9 ^ = cos 0 + 8+8 
The expansion and rearrangement of terms gives 


1 + ^ +0 {f 2 } 
R 

m -J 


8 = tan" 


cos | 

^ + 8 ) 

! 

AR 

1 + ir~ 

L m J 

- cos | 


s in J 

f 1 , 

— 0 + 8 

i 

1 1 

K; 

| + sin 

(M 


- tan" 


0 . s . 1 , ar n , s 

2 sin 8 sin — w cos — w + 8 

2 R l 2 

m \ 

2 cos 8 s in i i p + sin [ — 0 + 8 
2 r R 1 2 r 


o{f 2 } 


If 0 is not a small angle, 8 can be expressed as 


6 = tan 1 


-^ 5 - cosi 0-28 sin -0 + o{f 2 } 

R 2 ' 2 


2 sin — 0 [ 1 + + o { f 2 } 

2 l, . 2R 


or 


6 - - — 5 . co t i 0 - 8 + o{f 2 } 
2 R 2 


The use of equations A-23 and A-24 gives 


cos ( - 0 _ 8 


cos j - 0 + 8 + 8 


- 1 


A- 8 


(A-25) 


(A- 26) 


(A-27) 
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For large 0 when 9 = o { f }, 


c 

x 


cos — 0 + 9 sin-ii 

2 __ 2 

cos 0 - (9 + 8) sin — 0 


- 1 


= l+ 0tan-^i/> + (6' + S)tan — 0 + o{f 2 } -1 


= (29 + 8) tan i 0 + o{f 2 } 


(A-28) 


In conclusion, the use of these results in equation A- 21 gives 


where 


*a _ R s i n 0 

~2 ~2 


(1 + c x ) COS \' a 



and, for most practical purposes, 


6 - — ^5. cos— 0-8 + o{f 2 } 
2 R 2 


(A-29) 


(A-30) 


T he Horizontal Gradients of Potential 

The horizontal gradients of potential are obtained as follows. As V' is 
given by equation 33, 


V d = W 0 - U 0 " V a + ^ h d> 


A-9 
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differentiation with respect to x gives 


3x a 


Can 


3V a 

tan B 

3x a 




+ 7— tan/ 3 a 


(A-31) 


on expressing the differential as it appears in equation 39. As changes with 
respect to x relate to the spherop U = Up, it follows that (e.g., Heiskanen and 
Moritz, 1967, p. 313) 


where £ are the components of the surface deflection of the vertical in the 
directions x . 

a 

The horizontal gradient of normal gravity is zero in the direction of the x 2 
axis. The gradient along the x x axis is obtained from the formula for normal 
gravity (see entry against fi in section 1.2.1) as 

-^Z. - y /3 sin 2<f> /R = o{10 -5 mgal cm -1 }. 
e 1 c 


thus the term 


h, ^21 tan /S. -o{10 -1 tan /3 mgal} 

OXj 1 


and cannot be ignored in either mountainous regions or areas where the height 
anomalies have significant magnitude. 
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The term BV a /^x a is purely a function of the topographical gradient, the 
change being approximately 20 /j.gal for a 250 m change in elevation (IAG 1970, 
p. 72). This is three orders of magnitude smaller than the change in Ag and is 
not of significance in the context of the other terms involved. Thus for the 
present development 

By' -v 

- — tan /3 = - y £ tan fi + h _Z tan /3 (A-32) 

ox a <3x. 1 

a 1 

the second term being two orders of magnitude smaller than the first. 

B. Spherical Harmonic Expansions and the Geodetic Boundary Value Problem 
B.l The Product of Two Surface Harmonics 
Let 


and 


S = P (sin 0 ) (Ci cos m\ + C 2 sin m\) 

nm nm v 7 v A nm z nm 7 


S^ k = P ^ (sin 0 c ) (Ci 4 cos kA. + C 2 ^ sin k\) 


be two surface harmonics. The product 


S n m % =P nm ( sin ^ c ) P ^k( sin ^ c ) [C lnm C 4k C ° S!lA CO S k\ + C x nm C O S mA S i n k \ 

+ C 2nm C sin mA. cos kA. + C 2nm C 2 ^ sin mA. sin kA.] 


A-ll 



Ill 


= P„„ (sin 4> c ) P U (sin * c ) [I (C ln _ C[ u - C 2 ^ C'J cos (m + k) X 

+ 5 ( c Uu C2 »- + c ‘~- C U») sin (m + k) k + 1 ( c >- + c s». % „) cos ( "- k)X 

4 K C s.. - C '.. C ^) Sin <"- k)X ]' 


Further (e.g., Mather 1971c, p. 47), if sin <fi c - fj , , 


O) = — (1 -^ 2 ) 1/2ra 
2 n 


K 

Z> 


i) r 


(2 n - 2 r) 


r ! (n - m - 2r)!(n - r) 




n~ m - 2 r 


where k* = 1/2 (n - m) if (n - m) is even or 1/2 (n - m - 1) if (n - m) is odd. 
Hence 


P„„ p (! k M = 22 *■ ^ 

s = 0 

where p = n + i . It should therefore be possible to represent the above product 
by a relation of the form 


p i 

p„. oo w = Z Z Bi > Pii <rt ' 

i =0 j = 0 

It would suffice for the purposes of the present study to draw the conclusion that 
the product S nm S^ k can be fully represented by the set of surface harmonics 


A-12 



112 


CO 



p = 0 


where 


p 

S" = ^ P pq ( sin ^c) CA pq cos qX + B pq sin qk]. 

q= 0 

B.2 The Orthogonal Proper ty of Surface Harmonics over a Closed Surface 

It is well known that two surface harmonics of different degrees satisfy the 
relation 


S n d S = 0 if n / i. 

This orthogonal property is usually derived in the case when S is a sphere, the 
derivation being a consequence of the surface harmonic constituting a spherical 
harmonic term (e.g., Jeffreys and Jeffreys 1962, p. 636). This property could be 
extended to cover any continuous closed surface, e.g., the physical surface of the 
Earth. It must be assumed that the surface harmonic expansion exists and pro- 
vides unique definition at all points on the surface, the integration being taken 
over the element of solid angle between the appropriate limits defining the entire 
surface. This would imply that any set {0 c , k} would define a unique point on 
the surface. If sin 0 c = 

S =P ( a ) C, cos m \ + C„ sinm 

nm nm vr ^ ' 1 2 
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and 




= P 


k M [ c i 


X U 


cos k A. + Cj sin kA. 

£k J 


the surface integral 





cos m\ cos k A. + C, C' 
lnn > 


cos mA sinkA. 


+ C„ C' sin m \ cos k k + G, 
2 " m Uk 2 


C' 

2 4. 


sin m \ s i n k A. | d \x d A. 


= 0 if m^k, as 


/* 2rr 

r2rr 

/• 2tt 

1 cos m A. cos k A. d A. = 

I sin in A sink A d A = 

I cos raAsinkAdA = 0 

h 

Jo 

•'o 


If m = k, two terms remain, as 


/• 2tt 

I cos 2 m\dA = 
•'O 



1 + cos 2m\ 
2 


dk = 



1 - cos 2 m X 


d A. = TT . 


sin 2 m k d A. 


Therefore 




d a = 




Cm) d fj. . 
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The integral can be evaluated by using the relation (e.g., Hobson 1965, 
p. 99) 


p „„ oo = (- ')" 


(n + m) ! 
(n - m) ! 




O). 


and assuming that n> l , when 


f 

J -i 


P. M (M) <*- <-*<* + *) 

2 n+ ^ (n-m) ! n ! I ! 


f 1 d n ~' 

dfi n 


[(M 2 - l) n ] —- 7 — [<> 2 -lf] d/X 


d,u. 


/C + m 


On integrating the above equation (n - m) times by parts, 


f P n™ M 0*) d M (n + m - ) - - f (M 2 - l) n -1^-1 [(M 2 - l/] d/, 

J -i 2 n+,t n ! I ! (n - m) ! J - 1 df/' + a * 


0 as n > ■i, 


and the non- integrated product being zero on evaluation at the limits of integra- 
tion, due to always having a factor (m 2 - 1) (e.g., Mather 1971c, p. 51). 

If n<-t , then replace P^ m (/z) by ( _ m) (aO 311(1 proceed as before. 

It can therefore be concluded that in all circumstances of integration over 
a closed surface on which a set of surface harmonics provide unique definition, 



d ct = 


0 , 


n ^ % . 
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Thus the orthogonal properties of surface harmonics apply on integration 
over any closed surface which is single valued in the set {cp c , A.}. 

B.3 Laplace's Equation at the Surface of the Earth 

Attention has been drawn in note (i) to section 2.4, to the fact that the dis- 
turbing potential V d , as defined by equation 5, does not satisfy Laplace's equation, 
but the function V d , given by 


v* = v - V 

where V is the potential due to the Earth's atmosphere, does so at all points 
exterior to the Earth's surface. S in equation 29 is the physical surface of the 
Earth, in the strictest sense. 

The development in section 3 requires that V' be capable of representation 

d 

by a set of spherical harmonics at all points on S. This cannot be claimed to be 
the case if S is the exact physical surface of the Earth. On the other hand, V ' 
satisfies Laplace's equation at all points exterior to the Earth and right down 
to it. Thus V d can be expressed by a set of spherical harmonics at all points in 
space exterior to the Earth's surface and right down to but not on it, provided 
the reference ellipsoid is everywhere within the physical surface of the Earth. 
Another important point is that physical manifestations of V ' are not meas- 

d 

ured at the physical surface of the Earth but slightly exterior to it. Thus the 
surface S being defined in equation 29 is not the physical surface of the Earth 
itself, as a consequence of the nature of the observational data used to solve the 
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problem, the former being measured slightly exterior to the physical surface 
of the Earth. 

It can be concluded that physical validity exists for the representation of the 
scalar characteristics of the Earth's gravitational field by a spherical harmonic 
series), if based on the observations made "at" and exterior to the physical 
surface of the Earth. In the context of the present study which seeks resolution 
at the ±5 cm level, no significant errors occur if S is taken as a surface which 
is always slightly exterior to the physical surface of the Earth. Such a surface 
has the advantage that V' d satisfies Laplace's equation at all points on it and 
therefore physical validity exists for adopting the representation 

V * = r — — + o{f 2 }, n 1. 

L I Rn+l 

n = 0 

No common convention has been adopted for the definition of a "surface of 
measurement" for gravity determinations vis-a-vis the physical surface of the 
Earth. The accuracy requirements defined in section 4.3 for resolution of the 
geodetic boundary value problem at the ±5 cm level, call for the representation 
of the global gravity anomaly field so that systematic errors over large extents 
held to below ±50 /j. gal. This in turn requires the establishment of the global 
gravity control network with individual station accuracy at this same level. The 
surface of measurement, which is only relevant when defining the global control 
net, can be deviated from by individual observation stations by ±10 cm without 
affecting solutions of the boundary value problem to order e 3 . 
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C. A Non-Iterative Procedure for Evaluating the Gravitatio nal Terms in the 
Solution of the Geodetic Boundary Value Problem 

Section 4.2 shows that iterative techniques are necessary for solving the 
boundary value problem to order e 3 . This is not an economic procedure and it 
is compelling to search for a non-iterative solution. The basic boundary 
condition is 


V dp = V. 


ap 



X q tan lq 

r 2 




t an B 

~ a 



d a- (A-35) 


The first set of terms is obtained from equation A-29, while the second can 
be obtained from A- 9 as 


where 



(A-33) 


c x 3 - 8 sin <// - (1 + - c R ) cos \p + o {f 2 } (A- 3 4) 

The third term, defined by equation A-32, cannot be satisfactorily included 
in the major gravitational term, even though it has significant magnitude. The 
last term is given by equation 34, the quantity h rf from equations 11 and 25 as 

h d = 7 tV'+V, -(W 0 -U 0 )J (A-35) 
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and equation 43 can be written as 


where 


3 y 

Fh 


- — ( 1 + c ) 
R V 


(A-36) 


c^=f+m-3f sin 2 <fi c + o {f 2 } 


6 

(A-37) 


The combination of these relations gives 


V. = V 

dp a 


+ + 


(A-38) 


where 



- F (v/s h) 


3 h 


d cr 


(A-39) 


F - F (0, h) 



Sin ^( 1+C x )^- ( 1 + C x 3 ) 


(A- 40) 


and 



tan/3 a + h d 


t a 


n /3 1 


d cr 


(A- 41) 


On adopting a spherical harmonic representation for V the validity of which 
has been established in section B3, 



n i- 1 
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equation A- 3 9 can be written as 


Ia 2 


TT I r / untl 

J ^ n = 0 ' 


da n ^ 1 


(A- 42) 


The gravity anomaly is related to the disturbing potential through equations 

G 

13, A-35 and A-36 as 


Ag=- 


3 V * 2 V ' 


d h R 


C 1 + c 4>) + c A g 


(A- 43) 


where 


1 ,, Wo- U o V. 3 V, 

c A s = 2 8i t2 — — - 2 ~- 


R 


R Bh 


+ o {f 2 A g} (A- 44) 


Equation A-43 can be expressed in spherical harmonics as 


CO 



n = 0 


(n - 1 - 2 c ,) 


R 


n + 2 


+ C Ae 


n pi 1 


(A- 45) 


The spherical harmonic function as evaluated at the surface of the Earth, 
can be expressed a set of surface harmonics for the reasons given in section 3. 
Let 


n 1 = A g - (A- 46) 

If conventional practice is to be followed, it will be necessary to prove that 
the replacement 


Ag 


' = E g --E (n -‘- 2 

n=0 


c <i) 


R 


n + 1 


A-20 



120 


R n+ 1 G n 
n ~ 1 - 2c 0 


n =£ 1 


is valid, in which case I. can be written as 

w A 


(A-47) 


„ i rrR 2 v n + 1 

2 77 I I r Z_i n - 1 - 
J J n = 0 


+ F 
2 c 


G da, n / 1 


(A- 48) 


R 2 /r can be expressed as a set of surface harmonics for the reasons given 
in section B.l of the Appendix as 


00 



n = 0 


A non-iterative solution for the gravitational terms could only be obtained 
if I A can be transformed into an expression of the type 


I 


A 


K y-' n + 1 + F 

2ttZ-j 

n = 0 


S Ag do-, n^l. 

n c 


This would be possible only if F and c^ were unchanged on surface integra- 
tion. As this is not the case, it would appear that it is not possible to solve the 
geodetic boundary value problem without resorting to iteration when evaluating 
the expressions containing the gravitational terms. 
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D. The Error of Rep resentation of Gravity Anomalies at the Surface of the Earth 
The error of representation E{Ag> nm for an n° x m° area at the surface of 
the Earth is given by 


IX 

< e »so 2 = 2Z 

i=i 


(Ag. - Ag ) 2 

N 


(A- 49) 


where the Ag. are individual determinations of the gravity anomaly at N points 
within the area, and Ag is given by 


ag 4 E ag - 

; = n 


For a meaningful estimation, E {Ag) nm must be the mean of several such 
evaluations. Further, the gravity stations must be evenly distributed about the 
region center with N being very large. Estimates of the error of representation 
for various square sizes by several researchers are given in table A-l. Linear 
units have been converted to equivalent angular values using 10 km = 0.09 degrees. 
The figures given in table A-l are heavily, if not totally biased toward continental 
areas and with two exceptions, to regions where topographical gradients are 
small. It should be noted that no correlation is implied between the value of 
E {Ag} and elevation. Thus E {Ag} for an elevated plateau should have a mag- 
nitude similar to that for a coastal plain. In rugged mountainous terrain, E{Ag} 
can be 3 to 10 times as great, especially for smaller regions. This should not 
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however preclude the use of E{Ag} to represent the statistical characteristics 
of the global gravity anomaly field for error estimation purposes. For example, 
in such cases, E {Ag} 0 1 can be as much as 10 times greater than the value 
(±3 mgal) in flat areas, the effect being confined to 3% of the Earth's surface 
area where rugged topography occurs. This would increase e in equation 
116 by a factor of 10 while e NA in equation 117 will be twice as large. 
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Table A-l 

Error of Representation for Free Air Anomalies 


Source 

Square Size 
n = m 
(degrees) 

Latitude 

(degrees) 

E{A g} n 

mgal 

Region 


0 

0.0 



MO 

0.05 

0 

1.5 

USSR 

H 

0.1 

0 

2.8 

Finland 

MO 

0.1 

0 

2.8 

USSR 

H 

0.2 

300) 

5.4 

Global 

MO 

0.2 

0 

4.4 

USSR 

MO 

0.3 

0 

7.0 

USSR (Plains) 

MO 

0.3 

0 

10 

Urals 

MO 

0.3 

0 

25 

Caucasus Mountains, USSR 

H 

0.5 

300) 

9.0 

Global 

M 

0.5 

30* 

10.1 

Australia 

MO 

0.6 

0 

10.1 

USSR 

H 

1 

300) 

12.7 

Global 

H & M 

1 

45 

12.4 


M 

1 

30* 

13.5 

Australia 

MO 

1.1 

0 

13.3 

USSR 

MO 

1.6 

0 

16.0 

USSR 

H 

2 

300) 

17.6 

Global 

H&M 

2 

45 

20.8 


M 

2 

30* 

17.7 

Australia 

MO 

2.2 

0 

16.3 

USSR 

H 

5 

300) 

23.1 

Global 

H&M 

5 

45 

27.6 


H 

10 

300) 

26.6 

Global 

H&M 

10 

45 

29.3 
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Key: -1. Source H = Hirvonen 1956, p. 3. 

M = Mather 1967, p. 131 
H&M= Heiskanen & Moritz 1967, p. 279 
MO = Molodenskii et al 1962, p. 172 

Col. 3: (1) = based on global sample 

* = mean latitude for region of studies 
0 = converted from data for squares with equidistant sides 
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Table 1 

Range of 0 for Linear Variations in f(0) and F (0) to Order e 3 
N = Number of contributions t as. at equation 100 to the quadratures evaluation 
of N f 


Square size 
(in degrees) 
n 

f<0) 


F(0) 

Range of 0 
(in degrees) 

N 

(x 10 s ) 

Range of 0 
(in degrees) 

N 

(x 10 s ) 

0.001 

0 > 0.07 

10.0 

V 

o 

• 

o 

— 

0.005 

0 > 0.5 

0.6 

0 > 0.0 

— 

0.01 

0 > 0.8 

3.2 

o 

• 

o 

A 

-3- 

0.0 

0.05 

0 > 3 

0.4 

V 

o 

• 

CO 

0.0 

0.1 

0 > 6 

0.1 

V 

© 

• 

0.2 

0.2 

0 >13 

3.4 

0 > 2 

16.2 

0.5 

0 > 60 

2.3 

— 

— 


N t 


20.0 


16.4 
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Table 2 

Estimates of Systematic and Random Error Effects on the 
Computation of Stokes’ Integral 


n 

(in degrees) 

Maximum tolerable 
Systematic Error 
in A g Over Range n 
(± mgal) 

Error in N f due to 
E {Ag} n 
(± cm) 

0.01 

50 

0.03 

0.1 

5 

3 

0.5 

1 

50 

1.0 

0.5 

120 

5.0 

0.1 

1400 
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FIGURE CAPTIONS 
Figure 1. Gravity and Its Potential 

Figure 2. The Disturbing Potential at the Surface of the Earth 
Figure 3. The first Order Inertia Tensor of the Solid Earth 
Figure A-l. The Meridian Ellipse and the Topography 
Figure A-2. The Spherical and Ellipsoidal Coordinates 

Figure A-3. The x. Cartesian System in the Local Laplacian Triad and Geo- 
centric Spherical Coordinates 
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Figure A-2. The Spherical and Ellipsoidal Co-ordinates. 
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